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INTRODUCTION 

The work described herein was performed at the School of Aerospace 
Engineering, Georgia Institute of Technology during the period 12 February 
1986 - September 1988. Professors Erian A. Armanios and Lawrence W. Rehfield 
were the Principal Investigators. 

This research concerns the analysis and prediction of delamination damage 
that occur in composite structure on the on the sublaminate scale — that is 
the scale of individual plies or groups of plies. The objective have been to 
develop analytical models for mixed-mode delamination in composites. These 
includes: 

(1) the influence of residual thermal and moisture strains 

(2) local or transverse crack tip delamination originating at the tip of 
transverse matrix cracks 

(3) delamination in tapered composite under tensile loading. 

Computer codes based on the analytical models in (1) and (2) have been 
developed and comparisons of predictions with available experimental and 
analytical results in the literature have been performed. A simple analysis 
for item (3) has been developed and comparisons of predictions with finite 
element simulation is underway. 

The usual approach to dealing with localized phenomena is large scale 
numerical simulation and analysis, mostly by general purpose finite element 
codes. This approach is often supplemented by a "build and test" demonstra- 
tion, or series of demonstrations if repeated failures are encountered. While 
such approaches are often costly and inefficient, their major drawback is that 
fundamental principles are not discovered which provide the means to produce 
better results. Furthermore, the steps must be repeated all over again the 
next time a similar situation arises. 


Overview of the Research 


The research program can be separated into three elements: The influence 
of residual thermal and moisture stresses on the mixed-mode edge delamination 
of composites. The analyses of transverse crack-tip delamination and 
delamination analysis in tapered laminates under tensile loading. A detailed 
account of the analysis and applications of each element is provided in 
Appendices I through III. A brief description and summary of the major find- 
ings of each research element is presented in the following sections. 

Influence of Hygrothermal Stresses 

The sublaminate edge delamination analysis and code which had its origin 
in the research conducted under the earlier grant NAG-1-558 has been modified 
to include the effects of hygrothermal stresses. 

The model is applied to mixed mode edge delamination specimens made of 
T300/5208 graphite/epoxy material. Residual thermal and moisture stresses 
significantly influenced the strain energy release rate and interlaminar 
stresses. Both experienced large increases when thermal conditions were added 
to the mechanical strains. These effects were alleviated when moisture 
stresses were included. Thermal effects on the interlaminar shear stress and 
total energy release rate were totally alleviated for the same specific 
moisture content. Moreover, this value of moisture content was not signifi- 
cantly affected by the stacking sequence for the laminates considered. This 
work is presented in accomplishments 3,4 and 12. A complete derivation of the 
analytical model, Fortran program listing and applications are provided in a 
accomplishment 3 and Appendix I. 

Transverse Crack Tip Analysis 

Transverse crack tip delaminations originate at the tip of transverse 
matrix cracks. This situation appears in Figure 1 where a symmetric laminate 
made of 90° plies in the core region and angle plies in the top and bottom 
portions is subjected to a tensile loading. Under tensile loading transverse 


matrix cracks initiate in the core region reaching a saturation level at a 
crack spacing denoted by A in the figure. Delamination often initiate at the 
tip of these transverse cracks. This situation is depicted in the generic 
model shown in Figure 1 of a symmetric delamination growing from a transverse 
crack tip. 

Three analytical models, sublaminate shear, membrane and shear lag have 
been developed in order to estimate the saturated crack spacing distance. The 
saturation crack spacing corresponds to the distance from the crack where the 
broken plies regain their uniform stress/strain state i.e. where the 
interlaminar shear stress has decayed down to its far field (uniform) value. 
Based on the closed form expression for the interlaminar shear stress the 
crack spacing predicted by each model is presented in Table I. The experimen- 
tal result in the table is based on Reifsnider's work for a [0/90] s laminate. 
A complete derivation of these models is provided in Appendix II. 

The analysis of transverse crack tip delamination is presented in Appen- 
dix II and applied to [± 25/90 n ]s laminates in the range n=0.5 to 8 made of 
T300/934 graphite/epoxy material. Closed form expressions for the 
interlaminar stresses, total strain energy release rate and energy release 
rate components are obtained. A computer code based on this analysis is 
developed and implemented into an earlier mixed-mode edge delamination code 
developed under the previous NASA grant NAG-1-558 and presented in accomplish- 
ment 6 and 7. This code was used to estimate the critical strain levels and 
the associated delamination damage mode with increasing number of 90° plies in 
the [± 25/90 n ] s . Since mid-plane edge delamination is a possible damage mode 
in this type of laminates a mid-plane delamination analysis was developed and 
presented in accomplishment 10. A computer code based on this analysis is 
developed and implemented in the mixed-mode edge delamination code. The 
critical strain and associated delamination damage modes predicted appear in 
Figure 2 and Table II. The critical stresses and associated delamination 
damage mode are provided in Table III. 



Experimental results show that the local (crack tip) delamination phenom- 
enon is the predominant damage mode only for n=4, 6 and 8 specimens. For n<4 
edge delamination either in the mid-plane or in the 25/90 interface were 
observed in tests. The present analysis predicts mid-plane edge delamination 
for n=l/2 and 1 and mixed mode edge delamination for n=2 and 3, respectively. 
For n=4, 6 and 8 local delaminations are predicted to be the controlling 
damage mode with approximately 25 percent Mode II for the three specimens. 
The critical strains in Figure 2 and Table II are computed based on a fracture 
toughness values of 415 J/m , 140 J/m and 120 J/m for local delamination, 
mixed mode edge delamination and mid-plane edge delamination, respectively. A 
complete account of this work appears in Appendix II. 

Analysis of Tapered Composites 

A generic configuration of a tapered laminated composite is shown in 
Figure 3 where a 38 ply thick laminate is reduced to 26 ply by dropping three 
inner sets of plies. The basic analysis approach that is adopted utilizes two 
levels of modeling, a global scale and a local scale. The global scale is 
concerned with overall generalized forces and strains such as axial force and 
extension. A simple consistent deformation assumption is the foundation of 
this model. Global equilibrium equations are written and solved. 

The generalized strains determined from the global analysis serve to 
provide estimates for the key primary stresses in the belt of the tapered 
section. Local estimates of interlaminar stresses are determined on the basis 
of equilibrium condition. 

The total strain energy release rate is computed from the work done by 
the external applied loads. The work done by the external forces is based on 
the axial stiffness of the different elements in the tapered configurations. 
These elements are represented by the six sublaminates shown in Figure 4 where 
Ag denote the effective axial stiffness of the uncracked belt portion, the 
effective axial stiffness of the cracked belt portion. The uncracked belt 
portion in the tapered region makes an angle £ with the loading axis while the 



cracked portion makes an angle a due to delaminations along the taper and the 
uniform regions. These are denoted by a and b in Figure 4. The effective 
axial stiffness of the uncracked and cracked dropped plies are denoted by A u 
and A c respectively. The axial stiffness of the straight portion is denoted 
by A g for the belt and A^ for the core plies. 

The extent of delamination along the tapered and the uniform portion of 
the laminate has a significant influence on the axial stiffnesses A , A c and 
Ag^. This is due to the discrete number of ply drops in the core region as 
illustrated in Figure 5 and the pop-off of the delaminated belt region. 

A three-dimensional transformation is required in order to estimate the 
effective axial stiffness of the belt regions A g and A g ^. This is due to the 
belt layup and the orientation of the different belt portions to the loading 
axis as shown in Figure 6. 

The interlaminar stresses between the belt and the core plies are pre- 
dicted by considering the equilibrium of the belt region. The equilibrium 
equations are derived using a complementary potential energy formulation of 
the belt on an elastic foundation. The elastic foundation is made of two 
contributions. The first, is a continuous shear restraint provided by the 
resin pocket regions at the interface between the belt and the inner core 
plies. The second, is a discrete number of concentrated transverse normal 
(R.) and shear (T. ) forces at the ply drop locations as shown in figure 7 for 
i =1-4 . The distributed shear stiffness is denoted by G in Figure 8 while the 
transverse normal and shear stiffnesses at the ply drop locations are denoted 
by k^ and g. (i=l-4), respectively. 

The variation of the total strain energy release rate G with delamination 
a growing along the tapered region appears in Figure 9. The effect of 
delamination b along the uniform portion on a is also shown in the figure. 
The discrete jumps at a/h equal 20 and 40 correspond to the ply drop. A plot 
of the concentrated transverse normal and shear forces and the interface 
between the belt and the inner core appears in Figure 10. 



A detailed description of the analysis, closed form expressions for the 
total energy release rate and interlaminar stresses is provided in Appendix 
III. Additional refinements are planned within this general framework such as 
accounting for shear strains in the belt and increasing the number of 
sublaminate elements in the analysis. 
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Table I Comparison of Transverse Crack Spacing 


Model Saturated Crack Spacing 

(mm) 


Shear 

2 Sublaminates 

1.651 


4 sublaminates, a — * 0 

1.105 

Membrane 

1.004 

Shear Lag 

1.160 

Experimental 

1.131 















Table II Critical Strains and Associated Delamination Damage Modes 

Critical Strains (%) 

Number of Experimental Local Edge Delamination 

90° plies Delamination Mixed Mode Mid-Plane 

1/2 0.6058 1.6747 0.6819 0.6058 

1 0.5936 1.1685 0.6262 0.5677 

2 0.5934 0.8058 0.5964 0.6402 

3 0.5934 0.6427 0.5862 0.7582 

4 0.5369 0.5444 0.5810 0.8815 

6 0.3914 0.4264 0.5757 1.1133 


8 


0.3589 


0.3555 


0.5731 


1.3199 










Table III Critical Stresses and Associated Delamination Damage Modes 

Critical Stresses (MPa) 


Number of 

Experimental 

Local 

Edge Delamination 

90° plies 

Delamination 

Mixed Mode 

Mid-Plane 

1/2 

438 

1313.9 

535.0 

475.3 

1 

409 

784.0 

420.1 

3S0.9 

2 

324 

426.2 

315.4 

338.6 

3 

270 

285.1 

260.1 

336.4 

4 

211 

210.6 

224.7 

341.0 

6 

128 

134.7 

181.8 

351.6 

8 

94.4 

97.1 

156.6 

360.5 
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g. 1 Generic Crack-tip Delamination 
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Fig. 2 Critical Strains and Associated Delamination Damage Modes 
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Figure 3. Edge View of the Tapered Structure 



Figure 4. Modelling Approach 









Stiffnesses on Delamination 



Figure 6. Three Dimensional Coordinates in the belt Region 




Figure 7. Interlaminar Stresses 
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Figure 9. Variation of Total Strain Energy Release 
Rate with Delamination 
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Abstract 


Laminated composite structures exhibit a number of different 
failure modes. These include fiber-matrix debonding within individual 
layers, delamination or separation of the layers, transverse cracks 
through one or more layers and fiber fracture. These failures are 
influenced by enviromental conditions. Thermal and moisture 
conditions are significant factors in interlaminar delamination as 
well as the other modes of failures. 

A simple delamination analysis method is presented here. It is 
based on a shear -type deformation theory and includes hygro thermal 
effects. These enviromental conditions are applied to the strain 
energy release rate and interlaminar shear stresses. 

The method is applied to mixed mode edge delamination specimens 
made of T300/5208 graphite/epoxy material. Residual thermal and 
moisture stresses significantly influenced the strain energy release 
rate and interlaminar stresses. Both experienced large increases when 
thermal conditions were added to the mechanical strains. These 
effects were alleviated when moisture stresses were included. Thermal 
effects on the interlaminar shear stress and total energy release rate 
were totally alleviated for the same specific moisture content. 
Moreover, this value of moisture content was not significantly 
affected by the stacking sequence for the laminates considered. 



4 


Introduction 


Composites have been used in the aircraft industry since 1969. 
One aspect of concern for using composites in structures is separation 
of plies or delamination. This occurs in regions of stress raisers 
such as holes, ply terminations, cut-outs and free edges. 
Delamination along the free edge of laminates have been observed 
during testing and service. The presence of delamination, initiated 
by interlaminar stresses, causes redistribution of the stresses among 
the plies in a laminate. Thus, it usually results in a reduction of 
stiffness and strength. 

Figure 1 shows delamination in a rotor hub made of S2/SP250 

glass/epoxy. The specimen has delaminated in two places that can be 

depicted by the dark lines. Figure 2 shows delamination in a single 

cracked-lap-shear test specimen made of AS4/3502 graphite/epoxy [ 1 } . 

The specimen layup is [+45/0/90] quasi-isotropic with 8 plies in the 

6S 

strap and 40 plies in the lap. The tests were performed on a 
displacement controlled machine. Fiber glass tabs were attached to 
the specimen ends. Multiple, isolated free edge delaminations occur 
in the neighborhood of the lap/strap junction and the tab. Figure 3 
illustrates an I-beam section made of C3000/5225 graphite/epoxy woven 
cloth in the post-buckled regime. Free edge delamination depicted in 
the flanges precipitated the final failure in the specimen. Figure 4 
shows how delamination can take place in single cracked- lap - shear 
specimens subjected to compressive and tensile loading. Specimens A 
and B delaminated under compressive loading while C experienced a 
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tensile loading. These examples illustrate the importance of 
investigating delamination problems in composites. 

Thermal residual and moisture effects on a composite are 
practical enviromental conditions. Determining the response of these 
conditions on interlaminar stresses and energy release rates in 
laminated composites is the primary objective of this work. 

Delamination analysis can be based on two approaches. They are 
the strain energy release rate and interlaminar stresses. The 
interlaminar stresses are due to Poisson's ratio mismatch and 
difference in the coefficients of thermal and moisture expansion 
between plies. Delamination occurs when these stresses reach the 
interlaminar strength of the material. An alternative approach is 
based on the actual process of fracture rather than the strength 
concept. Delamination can propagate when the strain energy release 
rate at the crack front is sufficient to overcome the material's 
fracture resistance or toughness. 

The strain energy release rate can be obtained for three 
particular modes of crack action. These three modes are known as Mode 
I or opening mode, Mode II or forward shearing mode and Mode III or 
tearing mode. Several failure laws are formulated in terms of these 
modes { 2 ] . 

A simple analysis predicting interlaminar shear stresses and total 
energy release rate with the influence of thermal and moisture effects 
is developed. This simple approach is useful in understanding the 
basic mechanics of the problem and predicting the factors controlling 
the behavior. The method is directed to the needs of a practical 
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design environment. It is not intended to compete with large-scale 
numerical approaches, but rather to serve as the means for selecting 
and screening candidate configurations and providing trend 
information. Simple codes for a desktop computer have been created to 
analyze laminate configurations. 
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Literature Summary 

A historical discussion of previously developed work for 
predicting interlaminar stresses and energy release rates is presented 
to establish a basis for the proposed model and to permit the present 
work to be placed in proper perspective. 

Earlier analyses have reflected the prediction of interlaminar 
stress and energy release rate without hygrothermal conditions. 
0'Brien[3,4] investigated delamination onset and growth in 
graphite/epoxy laminates under uniform extension. A simple expression 
was developed for the total energy release rate based on classical 
lamination theory. Whitney and Knight [5] used classical laminated 
plate theory to develop an edge delamination specimen analysis. This 
work was limited to Mode I behavior. 

An analysis based on a shear deformation theory and a sublaminate 
formulation [6] was developed by Armanios and Rehfieldf 7 , 8 ] . This 
method provides good estimates for the interlaminar shear stresses. 
Energy release rate components are estimated based on these stresses. 
However, this method does not provide reliable estimates of peel 
stress since thickness strain is neglected. This analysis was limited 
to mechanical strain only. 

0 / Brien[9] modified his analysis to include thermal and moisture 
conditions. The influence of thermal effects was considered by 
Whitney ( 10] . 

The work of Reference 9 was based on a classical laminated plate 
theory. It was applied to mixed-mode edge delamination specimens. 
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The results were limited to strain energy release rate. Finite 
element modeling was used to determine the strain energy release rate 


energy release rate due to thermal effects. It decreased with the 

addition of moisture considerations . For a T300/5208 graphite/epoxy 

laminate with [±30/+3 0/9 0/93] layup, the thermal effect increased 

I s 

the total energy release rate by 170 percent when compared to 
mechanical loading alone. However, a specific moisture level of 0.75 
percent completely alleviated this increase. In calculating the total 
strain energy he showed that bending and coupling effects became 
important at high levels of moisture content. 

In Reference 10 a higher-order plate theory with transverse 
normal strain effects was developed. Peel as well as interlaminar 
shear stresses could be predicted by this method. The thermal 

influence on total energy release rate and interlaminar stresses was 
investigated using a Mode I specimen. Residual thermal effects showed 
a significant influence on the stresses and release rates. For a 
graphite/epoxy laminate of [O 3 , 9 O 3 ] layup , thermal effects increased 
the maximum peel stress by a factor of 2.7 over that of pure 


mechanical strains. 


In the present work both thermal and moisture influences are 
studied in a mixed-mode delamination specimen. The analysis includes 
total energy release rate as well as interlaminar stresses. 
Similarities between the interlaminar stesses and total energy release 
predictions with hygrothermal effects is investigated. 

In the subsequent sections, the analytical approach is developed. 
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The method is then applied to six graphite/epoxy laminates. A 
discussion of the hygrothermal effects on interlaminar shear stress 
and total energy release rate predictions is provided. 
Recommendations for further investigations are proposed. Appendices 
are included for completeness . The first provides detailed 
expressions of the governing equations. Appendix II defines the 
hygrothermal expressions and their use in the analysis. The last 
appendix shows a listing of the program used and sample output. 
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Analytical Approach 


Overview 

The sublaminate modeling approach describes the essential 
features of the laminate behavior in a simple way. A free edge 
delamination specimen is shown in Figure 5. A uniform strain , €, is 
applied in the axial direction. From symmetry, only one quarter of the 
specimen is considered. In Figure 6, the specimen is modeled as if it 
were composed of four distinct sublaminates. Sublaminates 2 and 3 
represent the group of plies above and below the crack, respectively 
in the cracked portion of the laminate, while sublarainates 1 and 0 
denote the same group of plies in the uncracked portion of the 
laminate . 

The use of sublaminates -- groups of plies that are conveniently 
treated as laminated units -- simplifies the analysis considerably. 
This approach is applied with confidence when the characteristic 
length of the response is large compared to the individual sublaminate 
thickness [6] . This sublaminate modeling approach has been verified in 
Reference 7 by comparison with a ply-by-ply finite element solution. 
These sublaminates are connected by enforcing the proper continuity 
conditions on stresses and displacements at their interfaces. 

Displacement fields within each sublaminate are defined as: 

u - x€ + U(y) + z(3 x (y) 

v - V(y) + z/? y (y) (1) 


w - W(y) 



11 


where u,v, and w denote the displacements relative to the x, y, and z 
axes, respectively. Shear deformation is recognized through the 
rotations JS X and £y. The governing equations for each sublaminate are 
derived using a virtual work approach. The derivation of the 
governing equations used in the development appears in Appendix I. The 
derivation is an extension of the work of Reference 8 with 

hygro thermal effects included. 

The constituitive relationships in terras of these force and moment 
resultants can be written as 

Nj-A^j +B ik »c k -N^ M (i,j,k - 1,2,6) 

Mj - Bjj € j + D jk < k - Mf (i,j,k - 1,2,6) (2) 

Qi - Ajj (i,j - 4,5) 

where the subscripts x, y, z, yz, xz, and xy are replaced by the 
subscripts 1-6 respectively. The force and moment resultants are 
denoted by N t and M ( - , respectively. Non-mechanical forces and moments 
resulting from the hygrothermal effects are labeled with a superscript 
NM. They are defined as: 

NM NM 

(Nj , M ) -/ Qij(ltZ) (Of. AT + b. C) dz (3) 

-h/2 J J 

The swelling coefficient is denoted by b j in Equation (3), the 
thermal coefficient by Oi- . The change in temperature is denoted bv 
AT and moisture weight gain by C. 
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The elastic stiffnesses A jj , By , and Dy are defined in terms of 
the sublaminate reduced stiffness Qy for a plane stress situation. 
These bear the classical definition. 

h/2 

< A ij • B ij • D ij > -jT /2 (l.z.z 2 ) Q j j dz (4) 

The equilibrium equations are: 

^xy,y + " t lx^ “ 0 

Ny,y + (t 2y - t ly ) - 0 


Qy,y + (P 2 

V 

- 0 



(5) 

**xy, y * Q x + 

h/2 

< c 2x 

+ t lx ) 

- 0 


M y,y ■ Qy + 

h/2 


+ t^ x ) 

- 0 


where t 2x , t 2 y, p 2 and 

c lx 1 

t iy* 

p^ denote the 

interlaminar s tress 


components in the x-z, y-z and z directions at the sublaminate upper 
and lower surfaces, respectively. These stress components appear in 
Figure 7. 

The displacements, resultant forces and moments, and interlaminar 
shear stresses in each sublarainate is governed by the displacement 
distribution (1), constituitive (2), and equilibrium (5) equations. 
These equations are applied to each sublaminate. The variables 


associated 

with 

each 

sublaminate are coupled through 

the 

continuity 

requirements 

at 

their 

interfaces . 

Enforcement of 

the 

boundary 

conditions 

lead 

to a 

solution for 

these variables. 

This 

procedure 



13 


discussed in general terms above is applied to the analysis of the 
edge delamination specimen shown in Figure 2 in the following sections. 

The response associated with sublaminates 1 and 0 shown in Figure 
2 is coupled through the continuity conditions at their common 
interface. The situation is different with sublaminate 2 and 3 where 
the continuity conditions are relaxed due to the presence of the 
crack. 

Uncracked Region: Sublaminates 0 and 1 

From symmetry conditions at the sublaminate bottom surface the 
shear stresses are zero. Interlaminar stresses at the top surface of 
sublaminate 0 are equal to those on the bottom of sublarainate 1 . 
Substituting these conditions into the equilibrium and constitutive 
relations and enforcing continuity of displacements at their common 
interface yields a homogeneous system of ordinary differential 

equations. These can be expressed in terms of the sublaminate 
rotations /J x and )3y. Assuming an exponential solution of the form 

(0ix*^ox »0iy >0oy ) - (0ix,£ox »0iy »^oy ) (6) 

results in a characteristic equation of the form 

EgS 8 + E 6 s 6 + E 4 s 4 + E 2 s 2 + E q - 0 (7) 

Parameter Eq depends only on the stiffness coefficients A^ 4 , A^ 
and A 45 for both sublaminates while Eg is predominantly influenced by 
the bending and coupling coefficients D - j and B,j . 


Thus, its 
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numerical value can be orders of magnitude smaller than the other 
coefficients. This results in the presence of a boundary zone in the 
response . 

The characteristic roots controlling the behavior of the laminate 
are determined from Equation 7 which has a closed- form solution. 

Crack Region of the Laminate: Sublaminates 2 and 3 

With this group of laminates, there are free surfaces at both the 
top and bottom of sublaminates 0 and 1 respectively. This is due to 
the presence of the crack. With the crack dividing the sub lamina te s , 
continuity conditions are not enforced at the boundary interface. 
This results in zero shear stresses at the surfaces of each 
sublaminate. Thus, the equilibrium and constitutive relations combine 
to produce a second order differential equation in terms of the 
sublaminate rotations #2X f° r sublaminate 2 and /? 3 y for sublaminate 3. 

'Interlaminar Stresses 

The arbitrary constants that are obtained from the eighth degree 
polynomial are determined by enforcing the stress free boundary 
conditions at the free edges of sublaminates 2 and 3, and the 
continuity of force, moment, displacement and rotations between 
sublaminates 0 and 3, as well as between 1 and 2. This yields the 
following expression for the interlaminar shear stresses. 

C X - "XYJ.Y- °j "' S,y (8) 

C Y - %,Y - T J 


( j “ 1 - 4 ) 


(9) 
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Parameters Tj and Gj represent the amplitude of the response. 

Energy Release Rate 

A complete formulation of the strain energy release rate appears 
in Appendix II. The total energy release rate can be determined by 
considering the work done by external forces. 

G 1/2 * dW/da (10) 

The total energy release rate is denoted by G and the crack 
length by a. The work done is defined as 

w - l/2 -C Nxi^mi d y < U > 

where subscript i denotes the sublaminate being referenced. The term 
represents the mechanical strain in each ply. This is defined as the 
difference between the total strain and the strain corresponding to 
free expansion for each ply. This strain is estimated by using a 
procedure similar to Reference 5. However, in Reference 5 the free 
expansion strain was determined by considering groups of plies in the 
cracked and uncracked regions of a Mode I edge delamination specimen. 
This approach is valid for a limited class of laminates. A 
ply-by-ply analysis rather than a sublaminate modeling should 
be used. In the following analysis, free expansion strains are 
determined on a ply-by-ply basis. 

From the symmetric conditions that exist in the uncracked section 
of the laminate, there exist no curvature. In the cracked portion. 
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the moment about the y-axis is assumed to be zero. Using both of 
these boundary conditions in Equation (2) yields the following. 

NM 

N Xl“ A ll € + A 12 € y + A 16 € z * N xl ^ ^ 

»,NM 

N xl“ A 11 €+ a 12 € Y + a 16 €z + B 12 K y" n x2 


Strain components € and G in Equation 12 are expressed in terms 

y z 

of the applied strain by 


„ ,NM 

Gy- C v € + C v 


(13) 


€ Z- C u €+ c 


NM 

u 


NM NM ^ . . 

The terms . C v , C u , C v and C u are functions of the extensional 

stiffness components of sublaminates 1 and 0. 


Using these expressions, Equation 12 can be re-written in the 

form. 


._k . k .^k _k . 

N xc - ( E c « c + T c ) 


(14) 


- ( E k € k + T k ) 
^XU ^ u X U * 


k k k k 

Parameters E c , T c , E u and T u are defined in Appendix II. 
Superscript k denotes the individual ply. Subscripts u and c 
represent the uncracked and cracked portions, respectively. The 
non-mechanical strain in each ply corresponding to a state of free 
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expansion is obtained by allowing the stress in each ply to vanish. 
This yields the following. 


k 

«u 


_k . 

- T u / 


E k u 


(15) 


k 


-Tc / Ec 


The strain corresponding to free expansion of the entire laminate 
is obtained by letting the resultant force vanish. The non-mechanical 
strain is 


€ NM - (t* - (tJ - T* ) 2a/b ]/| eS - ( eJ - E* ) 2a/b} (16) 


4c sic 4^ a k if k 

The terms Tj\ T c , E u and E* represent the summation of T u , T Q , E u 

k 

and E c over their respective sublaminates . These strain definitions 
for the effects of moisture and temperature can now be used in the 
general expression for the strain energy. The strain field is altered 
to represent the hygrothermal effects. The total strain for a 
sublaminate is expressed as: 


€ + € 


NM 


(17) 


The strain energy expression is given below showing the use of 
Equations (13) and (17). 

w - I ( € T + t£) (<= T - ) + ( Eu e T + t£ ) ( € T - ) j 


(18) 



Substituting this 
energy release rate 
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into Equation (11) yields the total strain 
per 


crack. 
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Results and Discussion 


Benchmark Study 


The analytical model is applied to the edge delamination specimen 
shown in Figure 5. The material considered is T300/5208 
graphite/epoxy . Its mechanical properties are listed in Table I. The 
coefficients of swelling and thermal expansions are also stated. The 
geometry of the specimen is given in Table II. Cure temperature for 
this material is 350° F. The ambient operating temperature is 70° F. 
The moisture level for all cases was allowed to vary from 0 to 1.2 
percent of the laminate weight. This moisture level reflects feasible 
conditions. The mechanical strain is taken as 0.00254. This 
particular value of strain was taken from test experiments [ 9 ] . It is 
considered a practical value for the material. 

Six laminates have been analyzed. They can be divided into two 
groups. The first group is composed of laminates [35/ -35/0^90^, 

I 35 /o / - 35 f 01 >- and t° / 35 / - 35 't 90 V la ” lna “ s are pr ° n<i co 

delaminate at the interfaces indicated by the arrow[9] . The Mode III 
in these laminates is negligible. This is due to the fact that 
relative sliding at the crack front and the interlaminar shear stress 
in the x-z direction, T Xz , is neglegible. The second group of 

laminates is [ 30/-60'75/15 L [ -35/55/10/-80]. and [ 35/80/- 55/- 10 ] . In 

f S' | s I s 

these laminates Mode I, Mode II, and Mode III are finite. The Mode 


III strain energy release rate component due to mechanical loading in 
these laminates are significantly large, ranging from 60 to 90 
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percent^ 7] . 


Interlaminar Shear Stresses 


The interlaminar stress at the delaraination interface appear 
in Figure 8-10, for the first group of laminates. The interlaminar 
shear stresses T XZ and Ty z for the second group of laminates appear in 
Figures 11 - 13. The labels M, M+T, and M+T+H stand for mechanical, 


mechanical and thermal, and mechanical, thermal and moisture 
respectively. 

The boundary layer of decay for all laminates ranged from 0.85 to 
0.93 of the laminate semi-width. In this context the boundary layer 
decay length is defined as the distance where the stress decays to 
5 percent of its maximum value. -The stress boundary zone is not 
significantly influenced by the environmental conditions. 

The magnitude of shear stress however showed a strong dependency 
on thermal and moisture conditions. At the delamination front, the 
ratio of stress with thermal effects as compared to pure mechanical 
loading ranged from 3.22 to 3.36 for the first group of laminates. 
This maximum was experienced at the crack tip. For the laminates 
where Mode III was present, this ratio ranged from 4.16 to 5.23 for 
Ty Z . The shear in the x-z direction showed a ratio of 1.4 to 2.16 for 
the maximum stresses. The maximum Ty Z stress for the second group of 
laminates was experienced at the crack front. However, the maximum 
T'xZ stress occured slightly ahead of the crack. This can be seen 
In Figures 11-13. 
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The addition of moisture alleviated the thermal effect. A 
moisture content of 0.4 has reduced the stress of thermal influence by 
approximatly 40 percent as compared to thermal influences alone. This 
trend is similar to the results of Reference 9. 

Numerical values of interlaminar shear stress at various moisture 
levels are provided in the sample output of Appendix III. 
Interlaminar shear stresses show numerical decrease with increase of 
moisture levels. 

Energy Release Rate 

The hygrothermal effect on total energy release rate appears in 
Figures 14 -15. The hygrothermal effects on total energy release rate 
show a similar trend to that of interlaminar shear stresses. Residual 
thermal stress tends to increase- total energy release rate while 
residual moisture alleviates this effect. The figures show that for 
total alleviation of the thermal effect, the specific moisture content 
ranges between 0.70 and 0.77 for all laminates. This indicates there 
exist a weak dependency on the stacking sequence. 

The effects of thermal conditions alone on the energy release 
rate does not correspond to the same numerical value as the 
interlaminar shear stresses. The total energy release rate of layups 
where Mode III was negligible showed a ratio of 5.1 for mechanical and 
thermal compared to mechanical conditions only. For the laminates 
where Mode III is finite, this ratio varied from 1.6 in the 

( 30/ - 60/75/- 15 ) layup to 3.37 in the [ 35/80/- 55/- 10 1 layup . 

T s I s 
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The total energy release rate in the first group of laminates is 
approximately the same for mechanical loading as shown in Figure 14. 
The influence of thermal and moisture does not appear to alter this 
trend. The energy release rate for the [ 35/- 35/j0/90] s laminate is 
indistinguishable from the [0/35^-35/90^ layup. The rate of 
alleviation due to moisture is the same for the three laminates. This 
is in contrast with the alleviation rate of the laminates where Mode 
III is finite as shown in Figure 15. For this class of laminates, the 
rate of alleviation due to moisture is different for each laminate. 

In some of the laminates, the rate of alleviation is not 
constant. There is a steep gradient in the rate of alleviation until 
the moisture content approaches the totally alleviated state. After 
such moisture content, the decrease in total energy release rate with 
respect to moisture addition is not as significant. 

■ It is worth noting there is a similarity between the strain 
energy release rate prediction and the interlaminar stresses for the 
totally alleviated state. This is shown in Figure 16 for a 
[ -35/55^10/-80] s layup. The specific moisture percent producing 
complete alleviation of the total energy release rate from the thermal 
effect is 0.76 as seen in Figure 15. The interlaminar shear stress 
distribution corresponding to this level of moisture is 
indistinguishable from the mechanical loading alone. The same 
conclusion was reached studying the other laminates. 



Conclusions 


A simple analysis was developed that predicted the influence of 
thermal and moisture effects on the interlaminar shear stresses and 
strain energy release rate. The analysis was applied to six 

mixed-mode edge delamination specimens. The results provide several 
significant findings. 

1. Residual thermal strain has a significant influence on the 
interlaminar shear stress and total energy release rate. 

The interlaminar stress and total energy release rate 
increased by 330 and 510 percent respectively over that of 
pure mechanical loading. 

2. Moisture tends to alleviate the thermal effect for both the 
interlaminar stress and energy release rate. At a specific 
moisture content of approximately 0.75 percent, the thermal 
influence is totally alleviated. 

3. The moisture content for total alleviation found from the 

total energy release rate analysis also produced an 

* 

interlaminar stress distribution similar to pure mechanical 
loading conditions. 

The first two findings are in agreement with the results of 
previous investigators. The third finding is new. It establishes a 
similarity in behavior between a delamination analysis expressed in 
terms of the energy release rate and the strength approach expressed 
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by the interlaminar stresses. 

These findings point to new directions for further inquiry. 
These are discussed in the following section. 
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Recommendations 


The thermal effects on the laminates showed a large increase in 
both the interlaminar shear stresses and strain energy release rate. 
The analysis should be supplemented with experimental tests to verify 
the result. Several fracture laws are expressed in terms of the 
strain energy release components, as well as the total strain energy. 
Further analysis should include predictions of these components in the 
presence of hygrothermal conditions. 

Throughout this work the temperature is assumed to be uniform 
through the thickness of the laminate. The same is true with the 
moisture. An approach corresponding to a practical environment method 
should account for temperature and moisture gradients in the laminate. 
In this situation, the hygrothermal gradients through the thickness 
may create an unbalance effect in an originally balanced construction. 
This consideration is of significant importance in aerospace 
structural components subjected to a large temperature difference 
between the upper and lower surface. 

The loading considered here is uniaxial extension. However, it 
is known that the load transfer points are not always in the plane of 
the laminate. Therefore, investigating laminate response under 
combined loads is' of great practical importance. It is recommended 
that bending, torsion as well as their combined effect be addressed. 

Findings by previous investigators suggested that delamination 
behavior in laminates subjected to fatigue loading follows static 
loading conditions. Further work is needed to investigate the 
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influence of hygrothermal conditions on the delamination of laminates 
under fatigue loading. 

Finally, the present analysis is applied to the mixed-mode edge 
delamination specimen. Extension of this work to other specimens such 
as the single- and double -crack- lap shear and the Mode II edge notch 
flexure specimen is recommended. 

When accomplished, these recommendations, together with the 
present research will provide a better understanding of the 
delamination problem in composites. Consequently, this will enable 
predicting, managing and ultimately preventing interlaminar fracture 
in laminated composites. 
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TABLE I - T300/ 5208 GRAPHITE/ EPOXY PROPERTIES 

£ 12 = 18.7 MS I 

E 22 = 1 .23 MSI 

G 12“ °* 8 3 2 MSI 

Poisson Ratio *= 0.292 

Swelling Coefficients of the Material direction: 
b ( 1 -di rection) ■ 0 
b ( 2 -direction) - 556 O b e/ %w eight 

Thermal Coefficients of the Material direction: 
tt(l-di rection) ■ -23 Ue/*f 
Qt(2-di recti on) • 1b«9uc/*F 


TABLE II - GEOMETRIC DIMENSIONS OF SPECIMEN 


Ply thickness « 0.005b inch 
Width = 1.5 inch 

Crack length ■ 6 x ply thickness = 0.032b inch 
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FIGURES 1-16 
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FIGURE 1 - GLASS /EPOXY ROTOR HUB AFTER DELAMINATION 
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FIGURE 2 - GRAPHITE/ EPOXY SINGLE CRACK- LAP- SHEAR SPECIMEN 
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FIGURE 5 - FREE EDGE DELAMINATION SPECIMEN 
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FIGURE 8 - COMPARISON OF INTERLAMINAR STRESS 
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FIGURE 10 - COMPARISON OF INTERLAMINAR STRESS 
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COMPARISON OF INTERLAMINAR STRESSES 



SHEAR STRESSES N KSI 


8 . 00 



FIGURE 12 


COMPARISON OF INTERLAMINAR STRESSES 


IN KSI 


10. 00 i 

o. oo - 
8 . 00 - 



FIGURE 13 - COMPARISON OF INTERLAMINAR STRESSE 
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FIGURE 14 - ENERGY RELEASE RATE DISTRIBUTION 
FOR LAMINATES WrTHOUT MODE ni 
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FIGURE 16 - TOTAL ALLEVIATED STATE STRESS DISTRIBUTION 


APPENDIX I 
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Derivation of the Governing Equations 

In this Appendix the governing equations for the sublaminate shown in 
Figure 7 are derived using the principle of virtual work. 

Consider a sublarainate of thickness h. The origin of a cartesian 
coordinate system Is located within the central plane (x-y) with the z-axis 
being normal to this plane. The material of each ply Is assumed to possess a 
plane of elastic symmetry parallel to xy as shown in Figure 6 . 

Stress and moment resultants are given below. 


(N X . Ny, N X y , 

,h/2 

Qx» Qy) " J^x. CT y» r xy, r x 2 » r yz) 

dz 


(H X # My • M X y) 

rV 2 

” 1 (^x* °y> f xy) z ^ 

■*h/2 


(i-D 

Because of 

the existence of a plane of 

elastic 

symmetry, the 


constitutive relations are given by 


*x 


Cll 


c x 

°y 


c 12 C 2 2 SYM 


c y 

°z 


c 13 c 23 c 33 


c z 

r xy 


c 16 c 26 C 36 C 66 


Yxy 


r rz 

. 

C 44 

SYM 


lyz 

r xz 


C 45 

c 55 


Yxz 


where C^j are components of the anisotropic stiffness matrix and 7 X y, 7y Z and 
7 xz are engineering shear strains. 

The displacements are assumed to be of the form 


u - U(x,y) + z0 x (x,y) 
v - V(y) + z0 y (x # Y) 
w - W(x,y) 


d-3) 
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where u t v and w are the displacement components in the x t y and z directions, 
respectively . Equation (1-3) in conjunction with the strain-displacement 
relations of classical theory of elasticity leads to the following kinematic 
relations 


c xx “ u *x* Z ^x,x 
€yy — V,y + z Py 9 y 

c zz “ 0 

7xy U*y + z (6 x ,y+ Py 9 x) 

7xz “ £x ♦ W,x 

7yz ~ 0y + w .y < X “ 4 ) 


Substitute Equation (1-4) into Equation (1-2) and put the results into 
Equation (1-1). This yields the following constitutive relations: 




A 11 

a 12 

a 16 

B n 

b 12 

N y 


a 12 

a 22 

a 26 

b 12 

b 22 

^xy 

— 

a 16 

a 26 

a 66 

b 16 

b 26 

«x 


B 11 

b 12 

b 16 

Dll 

d 12 

H y 


b 12 

b 12 

b 26 

Dll 

d 12 

M xy 


b 16 

b 26 

b 66 

Dl6 

d 26 

1 

r 


r 




Qy 


A 44 

a 45 

Qx 


a 45 

a 55 


where 

f h/2 

(Aij , By, Djj) - J C tJ (1, z, z 


b 16 


D.x 


N x 

b 26 


V.y 


Ny 

b 66 


U,r+V, x 

- 

N X y 

Dl6 


^x,x 


«x 

D26 


^y.y 


My 

d 66 


0x, + 0r,x 


«xy 


fiy + W t y 

0X + 


(1-5) 


and the non-mechanical terms are defined in Appendix II. 

The variation of the strain energy due to virtual displacements 6 u, 6v 


i 


and 6 w is 
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i 


r* 


i 


SV 


- [ (°X 


Sc x + CyScy + a z Sc z + r X y ^7 X y + Ty Z Sy yz + r xz 6y xz )&V 


d-6) 


where ic x , Scy , Sc z . 6 7xy» $7xz are t ^ ie strains associated with the virtual 
displacements. Using Equations (1-3) and (1-1) then integrating through the 
thickness gives 


SV 


- | [N x x + Ny5V,y + N X y 5U,y + Q x S0 X + Qy ( 60y + SW.y) + H X S 0 X t x 


+My 60 


y.y + M xy 


£ £x f y 1 dA 


d-7) 


The variation of the work done by the external forces and by the 
surface tractions is 


SW •» I (n x 5U + ny SV + q6V + m x 80 x + niy 80 y) dA 
J A 

+ f («; «U n + N ns 6V S + M n 60 n + M ns *0 s )ds (1-8) 

J S 

where a bar denotes values on the boundary. Variables n and s are coordinates 
normal and tangential to the edge , and 

n x ” c 2x * c lx 


n y - t 2 y - t ly 

q - P2 - Pi 

m x " | ( c 2x + c lx) 

m y “ “ (t2y + tly) (1-9) 


where n x and ny can be regarded as effective distributed axial forces. Terms 
m x and my are effective distributed moments and q is an effective lateral 
pressure. 


From the principle of virtual work the equations of equilibrium and 
boundary conditions are determined from the Euler equations and boundary 
conditions of the variational equation. 

SV - 6W 




(I-lCn 



51 


Substitution of Equations (1-7) and (1-8) into Equation (1-10) leads the 
following equations of equilibrium: 

N x,x + N xy,y + n x - 0 

**xy,x + Ny,y + Hy - 0 
Qx,x + Qy,y + Q “ 0 
**x,x + **xy,y - Q x + m x - 0 

^xy , x ^y , y “ Qy + ®y ”0 (1-11) 

and one member of the following five products must be prescribed on the 
sublaminate edges 


^n ^n » N ns u s * ^n^n » ^ns @s an< ^ Qn ^ 


( 1 * 12 ) 


For the ED specimen under uniform extension, U(x,y) in Equation (1-3) 
is given by 


U(x,y) - U*(y) + x€ 


(1-13) 


and the response is a function of y and z coordinates only. For this case the 
equilibrium equations (I -11) take the form 


N xy,y + n x - 0 

Ny ^ y ^ Hy ™ 0 

Qy,y + q - 0 
^xy,y ‘ Qx + m x “ 0 

My,y - Qy + my — 0 (1-14) 


Substitution of the constitutive relations in Equation (1-5) into Equation (I- 
14) yields the following equilibrium equations in terras of kinematic 
variables . 
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A 66*-yy 


U* 


n x 

A 26*-yy A 22^yy SYM 


V 


n y 

0 0 -A^Lyy 


V 

- - 

-q 

^66^yy ®26^yy ~ A 45*-yy (^66^yy * A 55) 


fix 


m x 

®26^yy ®22^yy “^ 44^-7 ^26^-yy " A 45 (^22^yy ~ A 44) 

t 

A 

1 

niy 


( 1 - 15 ) 

where the operators 


Lyy - d2/dy 2 

Ly - d/dy 

From these governing equations the basis of the work presented in this 
paper is formed. Appendix II gives a detailed formulation of the 
hygr c thermal terms and the formulation of the total strain energy release 
rate and interlaminar stresses. 
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APPENDIX II 
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Hvgro thermal Effects on E dee Delamlnation 


The displacement field and constitutive relations governing the free 
edge ply separation were presented in Appendix I. The hygro thermal 
expressions, represented with the superscript NM for non-mechanical, are 
defined as follows 


< Hf ) 



(II-l) 


where 

aj — Coefficient of thermal expansion 
— Swelling coefficient 
T — Local temperature 
T r — Reference temperature 
C - Specific moisture concentration 
Qij - Reduced stiffness coefficient 


The terms aj and are transformed as second order tensors with the 
assumption of no thermal or swelling shear strain. 

The concept of eublaminates is used when enforcing the boundary 
conditions. 


Cracked Sublaminates 
Sublaminate 2: 

The boundary conditions for this sublaminate are expressed as: 

Ny2 - N xy2 - My2 - Qy2 - 0 (II-2) 

M xy2,y * Qx2 “ 0 

Using the first three conditions in the governing equations, one can 
express V 2 , U 2 and /?2y * n terms of ^2x to obtain: 
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Sublaminate 3: 

The boundary conditions for this sublamlnate are given as: 
N y3 - N xy3 - 0 

M y3,y * Qy3 ** ® 

M xy3,y - Qx3 " 0 


(H-3) 


(II-4) 


These are used in a similar manner (as in sublaminate 2) to 




v 3,Y 

+ 

*26 


«3.Y 



A° R° R° 
x2 B 22 B 26 


.o _o _o 
A 16 B 26 B 66 


03, y 
03x,y 


N 


yo 


N 


xyo 


obtain 


- 0 


(II-5) 


These equations are then substituted back into the governing equations 
to obtain expressions for the force and moment resultants. They can be 
expressed in terms of the strain plus non-mechanical effects. 


UNCRACKED SUBUVHINATES 


Sublaminates 0 and 1: 


The boundary conditions of continuity at the interfaces must be 
satisfied. 


Ny Pb N xyl (0) - N y< £0>- N xyo <0)= 0 


(H-6) 


M X yl (0) - Mxy2(0) 
M yo (0) - H y3 (0) 
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" ^xy3(®) 


(11*7) 


M xyl (O) - M X y 2 (O) 

01x <0) “ ^2x(°> 

M y0 (O) - My 3 (0) (II*®) 

/?0y(O) ~ ^3y(0> 

M xy o(0) - M xy3 (0) 

^Ox(°) " ^3x(0> 

Enforcing equations (II-6) and (II-7) In the governing equations 
yields the following: 

f A 12 a 22 a 2s 1 f e 1 f"yil f"yu 1 


A 1 A 1 A 1 

A 12 A 22 A 26 


A 1 A 1 A 1 
16 26 66 


A° A ° A° 
A 12 A 22 26 


A° A ° A° 
16 A 26 66 


* N yU 

+ * N xylj s j G J “ 0 

N yU 

N xylj 


(II-9) 


fil 12 B 22 B 26 L 


A 1 B 26 B 66 


% n 


M yJj s j 

M xyij s j + B l s c a j 


(11-10) 


( j - 1- 4 ) 
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The expressions in (II-9) and (11-10) are defined below 
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( 11 - 14 ) 


Aj^ and A]_ are functions of Ay , By and Dy . The superscripts * 
implies a summation of the upper and lower sublaminates. 


Continuing with the derivation one can substitute the expressions set 
forth into equation (1-7) as well as 10 and 11 in the report. This gives the 
following expression for the total energy release rate. 
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( 11 - 15 ) 


The concept of free -expans ion in the x-direction is implemented to 
find the strain induced by the non-mechanical effects on the structure. 
Setting N x — O for each ply in Equation (1-5 ) and using the boundary 
conditions of (II-2), (II-4) and (II-6) allows the following. 
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where 

T k - h k cf q £ 2 ♦ h k cT - (C) 

E* - h k (Qn + Cv Q12 + Cu Qi 6 ) 

T k " Ql2 fl” h k + Q k 6 F^ h k -t- B k 2 F3** - (n2 M ) k 

E k - Q k x h k + Q k 2 h k Cd n + Q k 6 h k Cd k 2 + B 12 Cd 31 (II -17) 

NM 

Superscript k represents the ply. Expressions Cdy and are found 

by substituting the conditions of (II-2) into (II-3). 

r 1 A i „i V 1 f .1 -i 

a 22 a 26 b 22 a 12 b 26 

[cdl- - A26 a 66 b 26 a 16 b 66 

1 J 3x2 

b 22 b 26 d 22 b 12 d 26 

1 1 1 T* f 

a 22 a 26 b 22 N y 

f”” ’ “ a 26 a 66 d 22 ‘ N xyl 

j x 1 111 

b 22 b 26 d 22 l M yl 

Sublaminate 3 has By - 0 due to symmetry of the structure. When 

J NM NM 

considering these plies, the term F 3 and F2 are found by substituting the 

boundary conditions (II-4) into (II-5). 
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This gives a second expression of for this sublaminate 
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( 11 - 20 ) 


To find the total strain associated with the non-mechanical effects, 
it is necessary to sum the force over the entire structure and set it to zero. 
These are used in order to obtain Equation (16), (17) and (18) in the report 
on page 17. Substituting this in Equation (11-15) gives the total strain 
energy release rate expression per unit length 


G 


^.(^ T +6 (* T - 

2 K C 


*:> 


(B* 

u 


c T 4- I*) 
u 



) 


The expression c* - is 

c.u 

that ply. 


in essence the total mechanical strain of 


INTERLAMINAR STRESSES 

The interlaminar stresses of the structure are defined in Equations 
(8) and (9) 

fc x " ^xy,y ” ^xyij sj Gj e 

t y “ N y*y " N yU S J 2 G J e " sjy (j “ 3 ~ 4 ) (n-2i) 

While sj , the positive roots reulsting from the ploynominal 
Eg s 8 + Eg s$ + E 4 s 4 + E 2 s 2 + E 0 0, (11-22) 

are independent of the hygrothermal effects, the rest of the terms are not. 

Solving Equations (II-9) and (11-10) gives the term Gj while N X yjj an d 
Nylj are found from 
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(H-23) 


(j - 1,2, 3, 4) 


where vj , Uj and aj are found by imposing the boundary conditions on the mode 
shapes. They are dependent on the four values of sj as well as the 
sublaminate stiffness matrices. 



APPENDIX III 


PROGRAM START (INPUT, OUTPUT, TAPE5=«NPUT, TAPE 6 -OUTPUT) 

THIS PROGRAM IS FOR THE FINAL PAPER 8-16-87 
DIMENSION STATEMENTS 

REAL BG(C),E(9),GG(1|,C). MATR32(3,2), MATR33(3.3). MATR 3 ( 3 ) , 

C STRAIN (25), SAVE3 (3) » SAVE 33 (3, 3), WKAREA (99) , ZR 
COMPLEX SJ ( 8 ) 

DOUBLE PRECISION BNEG, A, C, DIFF1, DIFF2, UNSY (2) , UNSX (2) , 

C SSSS, SSSC, SSCC, SCCC, ZZZ0(0:50) ,J22,J26,FNM(3) .S1NM.S2NM, 

C MEMSY, MEMSX, F11M, F22M, SSI, SS 2 , SSY, SSX, ZZZ1(0:50), 

C TH I CK (1+0) .THETA (50) , El (50) . E2 (50) , CCCC, HSS (5) ,HSN1 ,HSN 2 , 
c Q ( 6 , 6 , 50 ) , ZO (0:40) , A 0 ( 6 , 6 ), A 1 ( 6 , 6 ) , U 13 , UU, BG 1 , BG 2 , 

C NXO(C), E15.E16.E17, E 18, E19, ZTT(0:C0), 

C NY 1 (C) , X (2) , Y (3) . CV, CU, W(2), ZZ (3) , VV 11 , VV12 
OOUBLE PRECISION ALPHA (C) , PHI (4), GAMMA (C) , NX1 (A) , 

C B 1 ( 6 , 6 ) , BO ( 6 , 6 ) , DO ( 6 , 6 ) , D1 ( 6 , 6 ) , F(C,C), VV13, WlC, J 66 , 

C RDLT , RTA1, RTA2 , RSB1, RSB2 , U12, Ull, A1NM, 

C NXY1(C), MY1(C), MXY1 (U) , WD (2,3) , CD (3, 2), WIDTH, 

C V12(50). V2 1 (50) , SS, CC, K 66 , K26, K22, 

C Z1 (0: AO) , FX, FY, G I (2,35) . K16.K12, H 66 

DOUBLE PRECISION SV(5), SU (5) . AL, SC, S ( 8 ) , DY, 

C G 1 2 (50) , G31 (50) , C2 , Cl. THETV, THETU, Gill (2,35). CS, 

C DEL, HO, HI, H22, HE, HG, HNY (50) , HNXY (50) , HM3.CVNM, 

C C11, C12, C22, C26, CCC, D, C55. C 66 , H 26 , SQ, DUM.CUNM, 

C CONY, CONXY, SMNY, SMNXY, 

C SB1 ,SB2,TA1 ,TA2,ATHMl (CO) ,ATHM2(C0) ,ATHM6(40) ,BSW2(C0) , 

C NMNYO , NMNXYO , G II (2 , 35) . 

C DVV11, DVV12, 0U13. OUlC, DF (C.C) , DX, ATH , CCONY , CCONXY , 

C BSW 6 (1*0) , DELTEMP, BSW1 (1*0) , CMOIST(35). 

C SIGX (0:A0, 79:120) , 

C NMNY1 , NMNXY 1 , NMMX 1 , NMMY 1 , S I G Y (0 : CO , 79 : 1 20) 

DOUBLE PRECISION NMSTO ( 50 ) , NMST2 (50) , NMST 3 ( 50 ) , TNC.UNCL, 

C T1 (50) , T12(50), T1 3 (50) , EX (50) , EXX (50) , EX3 (50) , JY, 

C EXNC, ESTAR, TSTAR, TNMST, GLC(0:50), NXNM(50) ,B12 (50) 

DATA Q/l800*0.0/, ZR/O.OO/ 

************************************************£************************ 

A************************************************************************ 

DATE OF PROGRAM : SEPT. 1, 1 987 

THIS IS THE FINAL PROGRAM FOR PREDICTING THE ENERGY RELEASE RATE OF 

COMPOSITE LAMINATES INCLUDING HYGROTHERMAL EFFECTS 

HOWEVER, IT ONLY CONSIDERS EXTENSION EFFECTS OF STRAIN WHEN DEALING 

WITH HYGRALTHERMAL EFFECTS LIKE WHITTNEY'S PAPER 

EXCEPT ON A PLY BY PLY ANALYSIS BASIS OF THE HYGRALTHERMAL EFF-ECTS 


THE INPUT ALLOWS FOR: THE LAMINATE LAY-UP TO CHANGE AND POSITION 
OF THE CRACK, DIFFERENT STRAIN VALUES TO BE EVALUATED, 
(UP TO AO LAMINATES AND 25 DIFFERENT STRAIN VALUES) 

AND FOR THE EVALUATION OF ONE MOISTURE CONSTANT OR A 
RANGE OF THE MOISTURE CONSTANT FROM 0 TO 1.2. 

THIS PROGRAM IS FOR THE GIVEN DATA TO BE IN ENGLISH UNITS. 

ALL LAMINATES ARE EVALUATED WITHOUT THRMAL EFFECTS AUTOMATICALLY 


************ ***************** ******************************************* 
*************************************** 3 ********************************* 


LZQ IS THE NUMBER OF DIFFERENT LAMINATE (OR CRACK POSITIONS) 
TO BE EVALUATED. 


RE AO (5,*) LZQ 

DO COO L ZZ = l.LZQ 
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OF POOR QUALITY 



RE AO (5 . *) WIDTH, NPLYO, NPLY1 , AL 
NEXTPL = NPLYO + 1. 

TPLY = NPLYO + NPLY1 

FOR EACH PLY IN THE SUBLAMINATE, THE MATERIAL CHARACTERISTICS 
MUST BE READ IN. 

PI = 1*. * ATAN(1.) 

HO - 0.0 
ATH *=0.0 
H 1=0.0 

DO 3 LK = 1, TPLY 
20 (LK) *= 0.0 
3 Z1 (LK) = 0.0 


DO 5 I - 1, NPLYO 

READ (5.*) THICK(I), THETA (I ) , El(l), E2(l) 

READ (5,*) VI 2(1) 

READ (5.*) G12(l), G3U0 
THETA ( I ) = THETA (I) * PI / 180. 

5 HO = HO + THICK (I ) 

DO 10 I « NEXTPL, TPLY 

READ (5.*) THICK(I), THETA (I) , E1(l), E2 (I) 

READ (5,*) VI 2 (I ) 

READ (5»*) G12 (I) , G3 1 (I) 

THETA (I) = THETA (I) * PI / 180. 

10 HI - HI + THICK(I) 

********************************************************************* 
********************************************************************* 
THESE ARE WRITE STATEMENTS TO CHECK THE INITIAL CONDITIONS OF THE 
SUBLAMINATE AND VALUES READ IN 

********************************************************************* 

********************************************************************* 

EACH PLY MAY HAVE DIFFERNT PROPERTIES SO THE PROPERTY OF EACH 

WRITE ( 6 , 289 ) 

CC 

WRITE (6,201) WIDTH 
WRITE (6,202) NPLY1, NPLYO 

WRITE (6,20*4) 

DO 15 J - 1, TPLY 
JJ = TPLY + 1 - J 
WRITE (6,206) J 

WRITE (6,207) THICK(JJ), THETA (J J) *1 80/P I 
WRITE (6,208) El (JJ) /1E+06 , E2(JJ)/lE+06 
WRITE (6,209) VI 2 (J J) 

15 WRITE (6,2114) G12 (JJ) /I E+06 , G31 (JJ) /1E+06 

******************************************************************** 
******************************************************************** 
DETERMINE THE 2 COMPONENT OF ALL LAMINATES 

CHECK = 0.00000001 
ZTT(O) =0.0 
Z0(0) = -HO / 2.0 
DO 20 I = 1, NPLYO 

ZTT(I) = THICK (I) + ZTT(l-l) 

20 20(1) = THtCK(l) + ZO(l-l) 


Z1 (NPLYO) = -HI / 2.0 
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DO 25 I = NEXTPL ,TPLY 

ZTT(I) - THICK(I) + ZTT(l-l) 

Z1 (I) - THICK(I) + Z1 (1-1) 

********************************************************************* 
FIRST READ IN THE NUMBER OF STRAINS TO BE EVALUATED AND THEIR VALUE 
THEN READ IN IF THE MOISTURE CONTENT SHOULD VARY OVER 0 TO 1.2 OR 
BE A CONSTANT. 

NSTRA = NUMBER OF VARIOUS STRAIN VALUES 

IF MOISTV - 1 ... CMOIST VARIES OVER 0 TO 1.2 
IF MOISTV = 0 ... CMO I ST IS A SPECIFIC VALUE 
********************************************************************* 
READ (5,*) NSTRA 
DO 27 J=l, NSTRA 
27 READ (5,*) STRAIN(J) 

DO 400 1ST - 1, NSTRA 
READ (5,*) RDLT.RSBl .RSB2.RTA1 ,RTA2 

WRITE (6,231) STRAIN (LST) ,RDLT, RSB1, RSB2, RTA1, RTA2 

READ (5,*) MOISTV 
IF (MOISTV. EQ.O) READ (5,*) CM 
IF (MOISTV. EQ.O) MMC - 1 
IF (MOISTV. EQ.O) WRITE (6,232) CM 
IF (MOISTV. EQ.1) MMC - 25 
IF (MOISTV. EQ.1) WRITE (6,233) 

********************************************************************* 
DO 300 JM - 1 ,MMC +1 

FIND Q'S AS WELL AS Q-BAR , SAVING Q-BAR 
AND READ AND CALCULATE THE HYGRO THERMO EFFECTS 

********************************************************************* 
DO 200 IZZ - 1,2 
III - 0 

IF (IZZ. EQ.1) JMM - JM 
IF (I ZZ.EQ.2) JMM - 0 
I F (JM.EQ. 1 . AND .IZZ. EQ.1) LIL - 1 
IF (JM.GT.l .AND. IZZ.EQ.2) GO TO 200 

NMNY1 - ZR 
NMNXY1 = ZR 
NMMX1 - ZR 
NMMY1 - ZR 

NMNYO - ZR 
NMNXYO « ZR 
HM3 - ZR 
SMNY - ZR 
SMNXY = ZR 

DO 24 1=1,5 
E(l) = ZR 
E (1+4) = ZR 
24 HSS(I) = ZR 

DO 26 1=1,6 
DO 26 J= 1,6 
NXNM ( I ) = ZR 
DF (I , J) = ZR 
AO(l,J) = ZR 
BO (I , J) = ZR 
DO ( I , J) = ZR 
A1 (I ,J) = ZR 
B1 (I , J) = ZR 
D1 (I ,J) = ZR 
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DO 28 MM * 1 ,TPLY 

28 I F ( THICK (MM) .GT. ATH ) ATH = THICK (MM) 

IZ2 = 2 IS FOR LAMINATE WITHOUT ANY HYGROTHERMAL EFFECTS 
IZZ * 1 IS FOR HYGROTHERMAL EFFECTS CONSIDERED 

DO 30 1=1, TPLY 

READ THE HYGROTHERMO EFFECTS, BOTTOM PLY IS FIRST AND UPWARD 

IF (I ZZ.EQ.2) GO TO 35 

IF (MOISTV.EQ.O) CMOIST(JM) = CM 
IF (M0ISTV.EQ.1) CMOIST(JM) = 0.05 * (JM-1) 
OELTEMP = RDLT 
SB1 = RSB1 
SB2 = RSB2 
TA1 = RTA 1 
TA2 ■= RTA2 
GO TO AO 

35 DELTEMP ■= ZR 

CMOIST(JM) = ZR 
SB1 ■= ZR 
SB2 = ZR 
TA1 = ZR 
TA2 = ZR 

AO V21 (I) - V12 (I) * E2 ( I ) / El (I) 

Cl 1 = El (I) / (1 - V12 (I) * V21 (I) ) 

C 1 2 = E2 ( I ) * VI 2 ( I ) / ( 1 - VI 2 (I ) * V2 1 ( I ) ) 

C22 = E2(l) / ( 1 - VI 2 ( I ) * V21 (I) ) 

CAA = G31 (I) 

C55 - G31 (I) 

C66 = G 1 2 ( I ) 

SS = DSIN (THETA (I )) * DS I N (THETA ( I ) ) 

CC = 1 - SS 

CS - 0.5 * DS I N (2*THETA (I ) ) 


SSSS = SS * SS 

SSSC = SS * CS 

SSCC = SS * CC 

SCCC = CC * CS 

CCCC = CC * CC 

Qd.l.l) - Cll * CCCC + 2 * (Cl 2 + 2 * C66 ) * SSCC 
c + C22 * SSSS 

QO.2,1) « (Cll + C22 - A * C66) * SSCC + C12 * ( SSSS 
C + CCCC ) 

Q (2, 2 , I ) = Cll * SSSS + 2 * ( C12 + 2 * C66) * SSCC 
C + C22 * CCCC 

Q(l,6,l) - (Cll - C12 - 2 * C66) * SCCC 
C + (C 1 2 - C22 + 2 * C66 ) * SSSC 

Q(2,6, I) = (Cll - C12 - 2 * C66 ) * SSSC 
C + (C12 - C22 + 2 * C66 ) * SCCC 

Q (6,6, I) « (Cll + C22 - 2 * C 1 2 - 2 * C66) * SSCC 
C + C66 * ( SSSS + CCCC ) 

Q(A,A, l)= CAA * CC+ C55 * SS 
Q (5.5. I) = CAA * SS + C55 * CC 
Q(A,5,I) = CS * (CAA - C55) 

Q (6.2. 1) = Q (2.6. I) 

Q(6,l,l) = Q (l ,6, l) 

Q (2, 1 , 1) = Q (1 ,2, I) 


HSS(l) = HSS(l) + Q (1 , 2, I ) 
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HSS (2) - HSS (2) + Q (2, 2 , I ) 
HSS (3) - HSS (3) + Q (2,6, 1) 
HSS (i*) « HSS(4) + Q (1 ,6, 1) 
HSS (5) « HSS (5) + Q (6.6, I ) 


ATHM1 (I) - TA1 * CC + TA2 * SS 
ATHM2 ( I ) = TA1 * SS + TA2 * CC 
ATHM6(I) * CS * ( TA2 - TA1) 

BSW1 (I) *= SB 1 * CC + SB2 * SS 

BSW2 ( I ) * SB1 * SS + SB2 * CC 

BSW6(I) « CS * ( SB2 - SB1 ) 

30 CONTINUE 

FIND THE A, B, AND D MATRICES FOR THE LOWER AND UPPER SUBLAMINATE. 
ALSO FINDS HYGRALTHERMAL EXPRESSIONS ON A PER LAMINA ANO PER 
SUBLAMINATE BASIS. 


ZZZO (0) = Z0(0) * Z0(0) * ZO (0) 
DO 45 I - 1, NPLYO 


NXNM(I) - (Q(l,1,l)*( ATHM1 (I) *DELTEMP + BSW1 (I) *CMO JST(JM) ) 

C + Q (1 ,2, I) * ( ATHM2 (I) *DELTEMP + BSW2 ( I ) *CMO I ST (JM) ) 

C + Q(l,6.l)*( ATHM6(l)*DELTEMP + BSW6 ( I ) *CMO I ST (JM) ) ) * THICK(I) 


NMNYO= NMNY0+(Q(1,2, l)*( ATHM1 ( I ) *DELTEMP + BSW1 (l ) *CMO I ST (JM) ) 

C + Q (2 , 2 , I) * ( ATHK2 ( I ) *DELTEMP + BSW2 (I) *CMOI ST (JM) ) 

C + Q(2.6, l)*( ATHM6 ( I ) *DELTEMP + BSW6 ( I ) *CMO 1 ST (JM) ) ) * THICK (I) 

NMNXYO* NMNXY0+(Q(1,6, l)*( ATHM1 (I ) *DELTEMP + BSW1 (I) *CMOIST (JM) ) 
C + Q (2,6, I) * ( ATHM2 ( I ) *DELTEMP + BSW2 ( I ) *CM0 1 ST (JM) ) 

C + Q (6,6, I) * ( ATHM6(I)*DELTEMP + BSW6 ( I ) *CM0 1 ST (JM) ) ) * THICK(I) 


ZZZO(I) = Z0(l) * ZO (I ) * Z0(l) 

B 1 2 (I ) = Q ( 1 , 2 , I ) *0 . 5* ( (ZO ( I ) *Z0 ( I ) ) - (ZO ( I - 1 ) *Z0 ( I - 1 ) ) ) 

DO 45 L -1,6 
DO 45 J *= 1,6 

I F ( REAL ( Q (J, L , I ) ).EQ.ZR) GO TO 45 
AO (J , L) - AO ( J , L) + Q (J , L, I ) * THICK(I) 

BO(J.L) - BO (J , L)+Q (J , L, I) *0.5* ( (ZO (I) *Z0 (I) ) - (ZO (I -1) *Z0 (I -1) ) ) 
DO (J , L) - D0(J,L)+Q(J,L. l)/3-0*( ZZZO(I) - ZZZO(i-l) ) 

45 CONTINUE 

ZZZ1 (NPLYO) « ZI (NPLYO) * Z1 (NPLYO) * Z1 (NPLYO) 


DO 50 I * NEXTPL.TPLY 

ZZZ1 (I) = Z1 (I) * Z1 (I) * Z1 (I) 


NXNM(I) = ( Q(1, !,!)*( ATHM1 ( I ) *DELTEMP + BSW1 ( I ) *CMO I $T (JM) ) 

C + Q(l,2,l)*( ATHM2 ( I ) *DELTEMP + BSW2 ( I ) *CMO I ST (JM) ) 

C + Q (1 ,6, 1 ) * ( ATHM6 (I ) *OELTEMP + BSW6 ( I ) *CMO I ST (JM) ))* THICK(l) 


NMNY1= NMNY1+ ( Q(l,2,l)*{ ATHM1 ( I ) *DELTEMP + BSW1 ( I ) *CMO I ST (JM) ) 
C + Q (2,2, I) * ( ATHM2 (I) *DELTEMP + BSW2 (I) *CMO I ST (JM) ) 

C + Q (2,6, I) * ( ATHM6(I)*0ELTEMP + BSW6 ( I ) *CMO I ST (JM) ) )* THICK(I) 



c 


NMNXY1- NMNXY1+(Q(1,6,I)*( ATHM1 (I ) *0ELTEMP + BSW1 (I ) *CMO I ST (JM) ) 
.C + Q(2,6,l)*( ATHM2 (I) *DELTEMP + BSW2 (I ) *CM0 1 ST (JM) ) 

C + Q(6,6,I)M ATHM6(I)*DELTEMP + BSW6 (I ) *CMO I ST (JM) ) )* THICK(I) 


NMMX1- NMMX1 + 0.5 * ( Z 1 ( I ) *Z 1 ( I ) - Z1 (I -1) *Z1 (I -1) ) * 

C ( Q (1.1, DM ATHM1 (I) *DELTEMP + BSW1 (I) *CM0I ST (JM) ) 

C + Q(1.2,l)*( ATHM2 (l)*DELTEMP + BSW2 (I ) *CM0 1 ST (JM) ) 

C + Q (1 ,6, I) * ( ATHM6(I)*DELTEMP + BSW6(I)*CM0IST(JM) ) ) 


NMMY 1 = NMMY1 + 0.5 * ( Z1 (l)*Z1 (I) - Z1 (I-1)*Z1 (1-1) ) * 

C ( Q(l,2,l)*( ATHM1 (I) *DELTEMP + BSW 1 ( I ) *CM0 1 ST ( JM) ) 

C + Q(2, 2, l)*( ATHM2 (I ) *DELTEMP + BSW2 (I ) *CM0l ST (JM) ) 

C + Q(2,6, DM ATHM6(I)*DELTEMP + BSW6 ( I ) *CM0 1 ST (JM) ) ) 


B 1 2 ( I ) - Q ( 1 , 2 , I ) *0 . 5* ( (Z 1 ( I ) *Z 1 ( I ) ) - (Z 1 ( I - 1 ) *Z 1 ( 1 - 1 ) ) ) 

DO 50 L«1,6 

DO 50 J«1,6 

IF ( REAl( Q(J,L, I) ).EQ.O ) GO TO 50 
A 1 (J, L) -A1(J,L) + Q(J,L,I) * THICK(I) 

B 1 (J , L) - B1 (J,L)+Q(J,L, 0*0.5* ((Zl (I)*Z1 (I))-(Z1 (1-1) *Z1 (1-1))) 
D 1 (J , L) - 01 (J,L)+Q(J,L,l)/3.0*( ZZZ1 (I) - ZZZl(l-l) ) 

50 CONTINUE 


****************** ************************************************** 


SEE IF COUPLING IS TAKING PLACE 

******************************************************** *********** 
COUPL = 2 
00 60 1-1.6 
DO 60 J*l,6 

IF ( REAL (BO (I ,J)) .GT. CHECK ) COUPL-1 
60 IF ( REAL (B 1 (I ,J)) .GT. CHECK ) COUPL-l 

IF (REAL (D 1 (2,6) ) .GT. CHECK) COUPL-1 
IF ( REAL (00(2,6)) .GT. CHECK ) COUPL-1 

IF ( COUPL. EQ.l .AND. LIL.EQ.1 ) WRITE (6,205) 

IF ( COUPL. EQ. 2 .AND. LIL.EQ.1 ) WR!TE(6,210) 
****************************************************************** 


CHECK THE SIGN OF THE PEEL STRESS *************************** 

HSN 1 - NMNY1 + NMNYO 
HSN2 = NMNXY1 + NMNXYO 

****************************************************************** 
HDD- HSS (2) * HSS(5) - HSS (3) * HSS (3) 

HE - HSS (3) * HSS (A) - HSS (1) * HSS (5) 

HE - HE /HDD 

HG - HSS (1) * HSS (3) - HSS (2) * HSS (A) 

HG - HG / HDD 
DO 65 l-l.TPLY 

HNY(I) = ATH * STRAIN(LST) * ( Q(l,2,l) + Q(2,2,l) * HE + 

C Q (2 ,6, I) * HG ) 

HNXY(I) = ATH * STRAIN(LST) * ( Q(l,6,l) + Q(2,6,l) * HE + 
C Q (6 , 6 , I ) * HG) 

65 CONTINUE 


DO 70 I = l.NPLYl 


UUUUNUU 


H«3 * HM3 + ATH * HNY(l) * ( NPLY1 - I + .5 ) 
SMNY « SMNY + HNY(I) 

70 SMNXY «= SMNXY + HNXY(I) 


IF ( HM3.GT.ZR) GO TO 85 

IF (LIL.EQ.1) WRITE(6,*)' CASE OF COMPRESSIVE PEEL STRESS 
WRITE (6,218) 

DO 75 l-l.TPLY 

WRITE (6,220)THETA(I) ,HNY(I) ,HNXY(I) 

WRITE (6,*) ' THE MOMENT CALCULATED WAS = ' ,HM3 

85 DO 80 1*1,6 
DO 80 J*1 ,6 

IF( ABS ( REAL (BO (I , J) ) ) .LT. CHECK) BO(l,J)«=ZR 
80 IF ( ABS ( REAL (B1 (I ,J)) ) .LT. CHECK) B1(I,J) -ZR 

******************************************************************** 
********************************************************* 

DEFINE SOME PARAMETERS NEEDED IN THE PROGRAM 

********************************************************* 

H22 - B1 (2.2) + HI / 2.00 * A1 (2,2) 

H26 - B1 (2,6) + HI / 2.00 * A1 (2,6) 

H66 - B1 (6,6) + HI / 2.00 * A1 (6,6) 

C22 - BO (2 , 2) + HO / 2.00 * A1(2,2) 

C26 - BO (2,6) + HO / 2.00 * A1 (2,6) 

C66 - 80(6,6) + HO / 2.00 * A1 (6,6) 

K22 * A 1 (2,2) + AO (2,2) 

K26 - A1 (2,6) + AO (2,6) 

K66 ■ A1 (6,6) + AO (6,6) 

K12 - A1 (1,2) + AO (1 , 2) 

K16 - A1 (1,6) + AO (1 , 6) 

D - K22 * K66 -( K26 * K26) 

E 15 - DO (2,2) - H0/2*B0 (2 , 2) 

E 1 6 - DO (2,6) - H0/2*B0 (2 , 6) 

£17 - DO (6,6) - H0/2*B0 (6,6) 

E 1 8 - BO (1,2) - H0/2*A0 (1,2) 

E19 - BO (1,6) - H0/2*A0 (1,6) 

VV11 - ( K26 * H26 - K66 * H22 ) / D + ( HI / 2.00 ) 

VV12 * ( K26 * C26 - K66 * C22 ) / D + ( HO / 2.00 ) 

VV13 - ( K26 * H66 - K66 * H 26 ) / 0 

VVU 4 « ( K26 * C66 - K66 * C26 ) / D 

U11 ■ ( K26 * H22 - K22 * H26 ) / D 

U1 2 *= ( K26 * C22 - K22 * C26 ) / D 

U13 *= ( K26 * H26 - K22 * H66 ) / D + ( HI / 2.00 ) 

U14 *= ( K26 * C26 - K22 * C66 ) / D + ( HO / 2.00 ) 

F (1,1) - 01(2,2) + B 1 (2,2) * HI / 2.0 + H22*VV11 + H26 * Ull 
F (2,1) * H22 * VV1 2 + H26 * U12 

F (3, 1) - D 1 (2,6) + B 1 (2,6) * HI / 2.0 + H22*VV13 + H26 * U1 3 
F (A, 1) - H22 * VV14 + H26 * U14 

F (2,2) = DO (2 , 2) - BO (2,2) * HO / 2.0 + C22 * VV12 + C26*U12 
F (3,2) = H26 * VV1 2 + H 66 * U12 

F (1+ , 3 ) = H26 * VV1L + H 66 * UU 

F (3,3) - D1 (6,6) + B 1 (6,6) * HI / 2.0 + H26*VV13 + H66 * U 1 3 

• F (4 , 2) = 00(2,6) - BO (2,6) * HO / 2.0 + C22*VV14 + C 26 * U14 
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c - 3 - 



ooo o ooo o ooo 


F(A,A) - DO ( 6 , 6 ) - BO( 6 , 6 ) * HO / 2.0 + C 26 ftVVlA + C 66 ft uiA 
DX - K22 * K 66 

DVV 11 - - K 66 * H22 / DX + ( HI / 2.00 ) 

DVV12 *= - K 66 * C22 / DX + ( HO / 2.00 ) 

DU13 « - K22 * H 66 /OX + ( HI / 2.00 ) 

DU 1 i* « - K22 * C 66 / DX + ( HO / 2.00 ) 

DF (1,1) *= D 1 (2,2) + B 1 (2,2) * HI / 2.0 + H22*DVV1 1 
DF (2,1) « H22 * DVV12 

DF (2,2) - DO (2 , 2) - 80(2,2) * HO / 2.0 + C22 * DVV12 
DF (A, 3) - H 66 * DU1I* 

OF(3,3) « Dl( 6 , 6 ) + B1 ( 6 , 6 ) A HI / 2.0 + H 66 * DU13 
DF (I* , i*) «= DO ( 6 , 6 ) - 80(6,6) * HO / 2.0 + C 66 * Dim 


W(1) « F(3,3)*( F (2,2) *F (A , I*) — F (A, 2 ) *F (4,2) )- F ( 3 , 2 ) *f ( 3 . 2 ) * 
C F(A,A) + 2*F(i»,3)*F(ii.2)*F(3.2) - F (2,2) *F (1* , 3) *F (fc, 3) 

W(2) « - F (3,3)*(F (2,2)*A0(5.5) + F (A,A) *A0 (I* , 1») - 2*F (I*. 2) 

* AO(l»,5) ) ~ A1 (5,5) * ( F (2,2) * F(A,A) - 
F (A, 2 ) * F (A, 2) ) + F(3,2) * F( 3 , 2 ) * AO (5.5) * 

2.0 *F(A,3)*F(3,2)*AO(A,5) + F(A,3) * F(A,3) * AO (A, A) 


X (1) -F (3,1) * (F (2,2) *F (A, A) -F (A, 2 ) *F (A, 2 ) ) -F (3,2) * (F ( 2 , 1 ) *F (A, A) 
~ F (A , 1 ) *F (A , 2) ) + F (A,3) * (F (2,1) *F (A, 2 ) - F (A, 1 ) *F (2, 2) ) 

X (2) « - F (3, 1) * ( F (2,2) * AO (5,5) + F (A, A) * AO (A, A) 

- 2 * F (A, 2 ) * AO (A, 5) ) - AT (A, 5) * ( F ( 2 , 2 ) * F(A.A) 

- F (A, 2) * F (A, 2) ) + F (3,2) * ( AO (5,5) * F (2,1) - 
F (A, 1) *A0 (A,5) ) - F (A , 3) * ( F (2,1)*AO(A,5) - F (A, 1) *A0 (A, A) ) 


Y(l) * F (3, 1 ) * (F ( 3 , 2 ) *F (A , A) - F(A, 3 )*F(A, 2 ))- F (3,3) * (F ( 2 , 1 ) * 
F (A, A) - F(A,1)*F(A,2))+F(A,3)*(F(2,1)*F(A,3)-F (A,1)*F (3,2)) 

Y (2) - 0 - A 1 (A,5)*(F ( 3 . 2 ) *F (A, A) - F (A, 2 ) *F (A, 3)) 

“ F(3,l) * ( F (3,2) *A 0 (5.5) " F(A,3) * AO (A, 5) ) + 

A1 (5.5) * ( F ( 2 , 1 ) * F (A, A) - F (A , 1) * F (A, 2) ) + 

F (3,3) * ( F (2 , 1 ) * AO (5,5) - F (A , 1) * AO (A, 5) ) 

Y(3) - A 1 (A,5) * (F (3,2) *A 0 (5,5) - F(A,3)*AO(A,5)) - A 1 (5,5) * 
c ( F ( 2 , 1 ) * AO (5,5) - F (A, 1 ) ft AO (A, 5 ) ) 


ZZ( 1 ) - F (3, 1) * (F ( 3 , 2 ) * F (A, 2 ) - F(A,3)*F(2,2)) - F (3.3) * (F (2, 1) * 

C F (A, 2) -F (A, 1 )*F (2,2)) +F ( 3 , 2 ) ft (F ( 2 , 1 ) *F (A,3)-F (A,1)*F (3,2)) 

ZZ ( 2 ) - F (3 , 1 ) * (F (A , 3) *A 0 (A , A) - F (3.2) *A0 (A,5) ) - A 1 (A , 5 ) * 

C (F (3.2) *F (A, 2 ) - F(A,3)*F(2.2))+ A1(5,5)*(F(2,D* F(A,2)- 

C F (A,1)*F ( 2 , 2 )) - F (3,3) * (F (A , 1 ) *A 0 (A , A) - F (2, 1) *A0 (A,5) ) 

ZZ (3) * 0 - A 1 (A , 5) * (F (A , 3 ) *A 0 (A , A) -F ( 3 , 2 ) *A 0 (A , 5 ) ) +A 1 (5,5) * 

C ( F (A, 1 ) ft AO (A, A) - F ( 2 , 1 ) * AO (A, 5) ) 

************************************************************************ 


NOW OBTAIN THE VALUES OF E SO THAT THE 8TH ORDER POLYNOMIAL MAY BE SOLVED 


C 

C 

C 


E (1) *= F (1,1) *W ( 1 ) - F (3.1)*X(1) + F (2, 1) *Y (1) - F(A,1)*ZZ(1) 

£ (3) = F ( 1 , 1) *W (2) - A 1 (A, A) ftW ( 1 ) - F (3, 1) *X (2) + A1(A,5)*X(1) 
+ F (2, 1) * Y (2) - F (A. 1) * ZZ (2) 

E(5) = (AO (A , A) *A0 (5,5) - AO(A,5)ftAO(A,5))*(F (1, 1)*F (3,3) 

- F ( 3 , 1 )*F ( 3 , 1 )) + (F (2 , 2 ) *A 0 (5,5) + F (A, A) *A0 (A, A) - 
2*F (A.2)*A0(A,5))ft(F (), 1)*A1 (5,5) - F (3. D *A 1 (A, 5) ) 



c - AT (It.lf) *W(2) + A1 (J.,5)*X(2)+ F (2,1)*Y(3) - F (4, 1) *ZZ (3) 

E(7) - - (AO (4,4)* AO (5.5) - AO (4.5) *AO (4.5) ) * (F (1 . 1) *A1 (5.5) 

C + F(3,3)*A1 (*.*)- 2*F (3. 1)*A1 (4,5)) " (A1 (4,4)*A1 (5.5)- 

C A1 (4,5) *A1 (4,5) ) * (F (2,2) * AO (5.5) + F (4, 4) * AO (4,4) 

C - 2 * AO (4,5) * F (4,2) ) 

E(9) = (AO (4,4) *AO (5,5) - AO(4.5)*AO(4,5) ) * 

C (A1 (4,4) * A1 (5,5) * A1 (4,5) * A1 (4,5) ) 

CALL UP SUBROUTINE TO SOLVE 8TH ORDER POLYNOMIAL 

NDEG - 8 
IER «= 0 

CALL ZPOLR (E.NOEG.SJ, IER) 

KK « 0 

******************************************************************** 

******************************************************************** 



IF 

(LIL.EQ. 

1) 

WR 1 TE (6 

,217) 



DO 

90 L * 

1. 

8 




S(L) - REAL 

(SJ (D) 




IF 

(REAL (SJ 

(L) 

|) .GT.O) 

KK - KK 

+ 1 

90 

IF 

(REAL (SJ 

(L) 

i) .GT.O) 

S (KK) - 

S(L) 


DO 95 kk 

SE 

1,4 



95 


IF (LIL. 

EQ 

.1) WRITE (6,221) 

KK, S (KK) 


******************************************************************** 


******************************************************************** 

NOW FIND THE UNCOUPLED S VALUES AND THOSE OF THE MEMBRANE 

BNEG- OF (1,1) * AO (4,4) + OF (2,2) * A1 (4,4) 

A « OF (1,1) * OF (2,2) - DF (2,1) * DF (2,1) 

C - AO (4,4) * A1 (4,4) 

SQ * DSQRT (BNEG * BNEG - ( 4.0 * A * C) ) 

DIFF1 - DABS ( BNEG - SQ ) 

DIFF2 = DABS ( BNEG + SQ ) 

IF (DIFF1 .GT.DIFF2) GO TO 100 

UNSY(l) = DSQRT ( (BNEG+SQ) / 2.0 / A ) 

GO TO 105 

100 UNSY(l) - DSQRT ( (BNEG-SQ) / 2.0 / A ) 

105 UNSY (2) = DSQRT ( (BNEG/A) - UNSY(l)* UNSY(1) ) 

BNEG* DF (3, 3) * AO (5.5) + DF(4,4) * A1 (5,5) 

A - DF (3.3) * DF (4,4) - DF(4,3) * DF(4,3) 

C - AO (5,5) * A 1 (5,5) 

SQ - DSQRT (BNEG * BNEG - ( 4.0 * A * C) ) 

0IFF1 « DABS ( BNEG - SQ ) 

DIFF2 « DABS ( BNEG + SQ ) 

IF (DIFF1.GT.DIFF2) GO TO 110 
UNSX(I) * DSQRT ( (BNEG+SQ) / 2.0 / A ) 

GO TO 115 

110 UNSX(l) = DSQRT ( (BNEG-SQ) / 2.0 / A ) 

115 UNSX (2) - DSQRT ( (BNEG/A) - UNSX(l) * UNSX (1) ) 

******************************************************************** 
IF (LIL.EQ.1) WRITE (6,224) UNSX(l), UNSX (2) 

IF (LIL.EQ.1) WRITE (6,223) UNSY(l), UNSY (2) 

******************************************************************** 


NOW THE S FOR THE MEMBRANE ONLY 



c 
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F11« • HI/2. * ( AT (2,2) * (0VV1 1 + DVV1 2*HO/H 1*A 1 (U , U) /AO (A ,U) ) 

+ B 1 (2,2) ) 

F22M - Hl/2. * ( A 1 (6,6) * (U 1 3 + UU»*H0/H1*A1 (5.5) /AO (5,5) ) 

C + B1 (6,6) ) 

MEMSY = DSQRT ( Al(l4,l*) / F 1 1M ) 

MEMSX «= DSQRT ( A1 (5,5) / F22M ) 

IF (LIL.EQ.1) WRITE (6,219) MEMSX, MEMSY 

******ft*******^******rt*****rtrt***>V*rt**ftftrt*rtrt*}V:fr*rtrt*ft*********rtrt*j'c*rt}'; 

OUM - AO (2,2) * AO (6,6) - AO (2,6) * AO (2,6) 

S1NM - ( AO(6,6)*NMNYO - AO (2,6) *NMNXYO ) / DUM 
S2NM - ( AO (2,2) *NMNXYO - AO (2,6) *NMNYO ) / DUM 
IF (LIL.EQ.1) WRITE(6,288) S1NM.S2NM 

IF ( COUPL.EQ.2) GO TO 130 

SOLVE FOR WD, CD, CU AND CV 

WD(l.l) - (AO (2 , 6) * AO (1 ,6) - AO (6,6) * AO 0.2) ) / DUM 

WD (1 , 2) - (AO (2,6) * BO (2,6) - AO (6,6) *B0 (2,2) ) / DUM 

WD(1,3) - (AO (2,6) *B0 (6 , 6) - AO (6,6) *B0 (2,6) ) / DUM 
WD (2, 1) * (AO (2,6) * AO (1 , 2) - A0(2,2) * AO ( 1 , 6) ) /DUM 

WD (2,2) - (AO (2,6) *B0 (2,2) - AO (2,2) * BO (2,6) ) / DUM 

WD (2,3) - (AO (2,6) * BO (2,6) - AO (2 , 2) *B0 (6 , 6) ) / DUM 


MATR33 (1,1) - A1 (2,2) 
MATR33 (1,2) ■= A1 (2,6) 
MATR33 0.3) - B1 (2,2) 

MATR33 (2,1) «= MATR33 (1,2) 
HATR33(2,2) = A1 (6,6) 
MATR33(2,3) - B1 (2,6) 

MATR33 (3.1) - MATR33 (1,3) 
MATR33 (3,2) - MATR33 (2,3) 
MATR33(3.3) * D 1(2,2) 


MATR32 (1,1) - - A1(1, 2) 

MATR32 (1,2) = - B1 (2,6) 

MATR32 (2,1) = - A1 (1,6) 

MATR32 (2 , 2) •= - B1 (6,6) 

MATR32 (3,1) = - B 1(1.2) 

MATR32 (3,2) - - D1 (2,6) 

IF (1ZZ.EQ.2) GO TO 122 

DO 120 1-1,3 
DO 120 K=l,3 
120 SAVE33 ( I ,K) = MATR33 0.K) 

SAVE3 (1) * NMNY1 
SAVE3 (2) = NMNXY1 
SAVE 3 (3) = NMMY1 
IRR = 0 

CALL LEQT2F (SAVE33. 1 . 3. 3. SAVE 3, 0, WKARE A, I RR) 
M=2 


122 
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N=3 
I A=3 
I RR=0 
I D0=0 

CALL LEQT2F (MATR33,M,N , I A, NATR32 , I DO .WKAREA , I RR) 
DO 125 1-1,3 

IF (I2Z-EQ.2) SAVE30) - ZR 
FNM(I) = SAVE3 (0 
DO 125 L=1 , 2 

125 CD (I , L) ~ MATR32(I,L) 


SC « DSQRT ( ( A1 (5,5) - Al(4,5) * A1 (4,5) / Al(4,4) ) 

C / ( D 1 (6,6) + B1 (2,6) *CD (1 ,2) + B1 (6,6) *CD (2,2) 

C + D1 (2,6) * CD (3.2) ) ) 

GO TO 135 

******************************************************************* 


IN CASE THE LAYERS ARE UNCOUPLED 


********************************************************************* 


130 


c 


DY - -1/ ( A1 (2,2)*A1 (6,6) - A1 (2,6) * A1 (2,6) ) 

CD (1 , 1) - ( A1 (6,6)*A1 (1,2) - A1 (2,6)*A1 (1,6) ) / OY 

CD (2, 1) - ( A1 (2,2)*A1 (1,6) - A1 (2,6) * Al(l,2) ) / OY 

CD (3 , 2) - ( A1 (2,2) *A1 (6,6) + A1 (2,6) *A1 (2,6) ) * 

D 1 (2,6) / D 1 (2,2)/ DY 
CD (1,2) - ZR 

CD (2,2) - ZR 

CD (3, 1) = ZR 


DR - -1 / ( AO (2,2) *A0 (6,6) - AO (2 , 6) *A0 (2 . 6) ) 

WD (1,1) - ( AO(6,6)*AO(l,2) - AO (2 , 6) *A0 (1.6) ) / OR 
WD(2,1) - ( AO (2, 2) *A0 (1,6) - AO (2,6) *A0 (1 ,2) ) / DR 
WD (1 , 2) - ZR 
WD(1,3) “ ZR 
WD (2 ,”2) - ZR 
WD (2 , 3) - ZR 

SSY - DSQRT ( AO ( 5 , 5 ) / DO (6,6) ) 

SSX - DSQRT ( AO (4,4) / DO (2,2) ) 


FNM(1) - (A1 (2,6) *NMNXY1 - A1 (6,6) *NMNY1) / OY 
FNH (2) - (A 1 (2,6) *NMNY1 - A1 (2,2) *NMNXY1 ) / DY 
FNM(3)«(A1 (2,6)*A1 (2,6) - A1 (2,2) *A1 (6,6) ) /DY * NMY1/D1 (2,2) 

SC « DSQRT ( (A1 (5,5)* A1 (6,6)) / (A 1 (6 , 6) *D 1 (6, 6) 

C - (B1 (6,6) *B1 (6.6) ) ) ) 


************************************************************** 

****************************************************************** 


135 Cl « B1(l, 6) + CD (1 , 1) *B1 (2,6) + CD (2, 1) *B1 (6,6) 

C + CD (3. 1) * D 1 (2,6) 

C2 = D1 (6,6) + CO (1,2) *B1 (2.6) + CD (2,2) *B1 (6,6) + 
C CD (3,2) * D 1 (2,6) 


J22 = 00 (2,2) + BO (2,2) * WD(1,2) + BO (2,6) * WO (2,2) 

J66 = DO (6,6) + BO (2 , 6) * WD(1,3) + BO(6,6) * W0(2,3) 

J 26 = DO (2 , 6) + BO (2 , 2) * WD(1,3) + BO(2,6) * WD(2,3) 

BNEG = J22 * AO (5,5) + J66 * AO (4,4) - 2. * J26 * AO (4,5) 

A = J22 * J66 - J26 * J26 



C - AO (4,4) * AO (5,5) - AO (4,5) * AO (14,5) 


SQ •= DSQRT ( (BNEG * BNEG) - I 4 .O * C * A ) 

OIFF1 - DABS ( BNEG + SQ ) 

DIFF2 *= DABS ( BNEG - SQ ) 

IF (DIFF1 .GT.DIFF2) GO TO 140 
SSI = DSQRT ((BNEG + SQ) / 2. / A ) 

GO TO 145 

140 SSI -= DSQRT ( (BNEG-SQ) / 2. / A ) 

145 SS2 « DSQRT ( (BNEG/A) - SSI * SSI ) 

SSY «= DSQRT ( AO (4,4) * J22 ) 

SSX ■= DSQRT ( AO (5,5) * J66 ) 

******************************************************************** 


IF (LIL.EQ.1) WRITE (6,*) ' SI AND S2 * ' , SS 1 , SS2 

IF (LIL.EQ.1) WRITE (6,*)' SX AND SY = \SSX, SSY 

******************************************************************** 


CVNM - ( K66 * (NMNY1 + NMNYO) - K26 * (NMNXY1 + NMNXYO) ) / D 
CUNM - ( K22 * (NMNXY1 + NMNXYO) - K26 * (NMNY1 + NMNYO) ) / D 

CV - STRAIN (LST) / D * ( K26 * Kl6 - K66 * K12 ) + CVNM 
CU «= STRAIN (LST) / D * ( K26 * K12 - K22 * Kl6 ) + CUNM 

****************************************************************** 


NOW FIND SOME OF THE NEEDED CONSTANTS 

FIRST DO LOOP IS TO VARY THE VALUES OF S 


DO 150 I = 1,4 

FORM THE A MATRIX (MATR32) 

MATR33 (1,1) = - ( 
MATR33 (1,2) = - ( 
MATR33(1,3) - ~ ( 
MATR33(2,1) ■= - ( 

MATR33 (2, 2) = - ( 

MATR33 (2,3) - “ ( 

MATR33 (3. 1) - ~ ( 
MATR33 (3,2) - ‘ - ( 
MATR33 (3,3) - - ( 

MATR3 (1) - (F (1,1) 
MATR3 (2) - F (3, 1) 
MATR3(3) - F (2,1) 


AND ITS B (MATR3) 


F ( 3 , 1 ) 

A 

SO) 

* 

SO) 

- 

A 1 (4,5) 

) 

F (2,1) 

* 

SO) 

A 

SO) 

) 



F (4,1) 

* 

s CD 

A 

SO) 

) 



F (3,3) 

* 

SO) 

A 

s(i) 

- 

A l (5,5) 

) 

F (3,2) 

* 

S(l) 

A 

SO) 

) 



F (4,3) 

A 

S(I) 

A 

so) 

) 



F (3.2) 

A 

S(I) 

A 

so) 

) 



F (2,2) 

A 

s(i) 

A 

SO) 

- 

AO (4,4) 

) 

F (4,2) 

* 

SO) 

A 

SO) 

- 

AO (4,5) 

) 

* S (1) 

A 

s(i) 

) 

- A 1 (4,4) 


S(l) 

* S(l) 

- 

A 1 (4,5) 




S ( I ) * S ( I ) 


CALL UP ROUTINE TO FIND THE VALUES OF ALPHA, PHI AND GAMMA 


M=1 

N=3 
I RR=0 
I 00=0 

1A-3 

CALL LEQT2F (MATR33«M,N, I A.MATR3, I DO . WKARE A , I RR) 

ALPHA (I) = MATR3 (1) 

PHI (I) = MATR3 (2) 

GAMMA (I) = MATR3 (3) 

SV(I) = VV11 + ALPHA (l)*VV13 + PHI (l)*VV12 + GAMMA ( I ) *VVl4 
150 SU ( I ) = Ull + ALPHA ( 1 ) *U 1 3 + PHI (l)*U12 + GAMMA(l)*Ul4 


00 155 I - l,A 


155 


NX1 (I)«A1 (1,2)*SV(I) + A1 (1,6)*SU(I) + 
NY1 (l)-Al (2,2)*SV(I) + A1 (2.6) *SU (I) + 
NXY1 (I)*=A1 (2,6) *SV(1) + A1 (6,6)*SU(I) + 
MY1 (I)=B1 (2,2) *SV(I) + B1 (2,6)*SU(I) + 
MXY1 (I)«B1 (2,6) *SV(I) + B1 (6,6) *SU(I) + 
NXO(I) ®A0 (1 ,2) *SV (l) + A0(1,6)*SU(I) + 
GG (1 , I ) = NY 1 (I) 

GG (2, 1) « NXY 1 (I) 

GG (3, 1 ) - MY 1 (I) 

FTH * C2 * SC 

GG (A, I) *= MXY1 (I) + FTH * ALPHA ( 


B1 (1,2) + B1 (1 , 6) *ALPHA ( I ) 
B1 (2,2) + B1 (2,6) *AIPHA (I) 
B1 (2,6) + B1 (6,6) *ALPHA (I) 
D1 (2,2) + D 1 (2,6) *ALPHA (I) 
01 (2,6) + 01 (6,6) *ALPHA (I) 
E 18*PH I (I) + E 1 9*GAMMA ( I ) 


)/ s CD 


A 1 NM - B 1 (2,6)* FNM(l) + B1 (6,6) *FNM (2) + 01(2,6) * FNM (3) 


BG (1) - A1 (1,2) *STRAIN(LST) + CV * A1 (2,2) + A1 (2,6) * CU-NMNY1 

BG (2) « A 1 (1,6) *STRA I N (LST) + CV * A1 (2,6) + A1 (6,6) * CU-NMNXY1 

BG (3) - B1 (1 ,2) ^STRAIN (LST) + B1 (2,2) * CV + B1 (2,6) * CU-NMMY1 

BG (A) - B1 (1,6) *STRA I N (LST) + B1 (2,6) * CV + 61 (6,6) * CU-A1NM 

C - Cl * STRAIN (LST) 

BG 1= BG (1) 

BG2 - BG (2) 


M=1 

n=a 

I A=A 
I DD=0 
1 RR=0 

CALL LEQT2F (GG ,M, N, I A, BG , I 00 .WKAREA, 1 RR) 
*************>v***********5ftr***yc*********^*****jv*Ajv***/f**************yt 

TVNM = S1NM - FNM(l) + HI / 2.0 * FNM (3) 

TUNM = S2NM - FNM (2) 

THETV * -CD (1,1) + WD (1 , 1) + H 1/2 .0*CD (3, 1) 

THETU - -CO (2,1) + WD (2 , 1 ) 

THETV = THETV + TVNM 
THETU ■= THETU + TUNM 

IF (LIL.EQ.1) WRITE (6,215) THETV, THETU 
IF (LIL.EQ.1) WRITE (6,216) SMNY, SMNXY 

THE STEPS USED TO FIND THE TOTAL ENERGY RELEASE FROM USING A PURE 
EXTENSION ANALYSIS FOR THE HYGRALTHERMAL EFFECTS. (SIMILAR TO 
WHITNEY'S). ANALYSIS 1$ CARRIED OUT ON A PLY BY PLY BASIS 


ZV - ( K26 * K 1 6 - K66 * K12 ) / 0 
ZU = ( K26 * K12 - K22 * Kl6 ) /D 
C 

DO 162 LL ■= l.TPLY 
C 

c 

EX (LL) = THICK(LL) * ( Q(1,1,LL) + ZV*Q(1,2,LL) + 

C ' ZU*Q(1,6,LL) ) 

Tl(LL) = NXNM(LL) - CVNM*TH I CK (LL) *Q (1 , 2, LL) - 
CUNM * THICK (LL) * Q(1 ,6,LL) 

EXX (LL) = Q(1,1,LL)*THICK(LL) + Q (1 , 2 , LL) *TH I CK (LL) *CD (1 , 1) 


C 



c 


+ Q(1,6,LL)*THICK(LL)*CD(2,1) + B12(LL)*CO(3,l) 

T1 2 (LL) ■= NXNM(LL) - FNM (1) *Q (1 , 2, LL) *TH I CK (LL) - 

C FNM (2) *Q (1 ,6, LL) *TH I CK (LL) - FNM (3) *B 1 2 (LL) 

EX3(LL) = Q(1,1,LL)*THICK(LL) + WO (1 , 1) *Q (1 ,2, LL) *THI CK (LL) 
C + WD(2,1)*Q(1,6,LL)*THICK(LL) 

T1 3 (LL) = NXNM(LL) - Q (1 , 2 , LL) *TH I CK (LL) *S1NM - Q(1,6,LL) 

C * TH I CK (LL) * S2NM 

IF (I ZZ.EQ.2) NMSTO(LL) « 0.0 
IF (I ZZ.EQ.2) NMST2 (LL) = 0.0 

IF (I ZZ.EQ.2) NMST3 (LL) = 0.0 

IF (I ZZ.EQ.2) GO TO 162 


NMSTO(LL) = T1 (LL) / EX (LL) 

NMST2 (LL) = T1 2 (LL) / EXX (LL) 

NMST3 (LL) = T1 3 (LL) / EX3(LL) 

162 CONTINUE 

WRITE (6,*)' JMM, NMSTO, 2. 3 OF ALL PLYS ',( NMSTO (JP) , 
C NMST2 (JP) ,NMST3(JP) , ' — ' ,JP«1,TPLY) 

WRITE (6,*)’ EX, EXX EX3 OF ALL PLYS *,( EX (JP) , 

C EXX(JP) ,EX3(JP) 1 , JP=1 ,TPLY) 


163 


TNC - 0.0 
EXNC = 0.0 
TSTAR =0.0 
ESTAR = 0.0 
DO 163 LK = l.TPLY 
TNC = TNC + T1 (LK) 

EXNC = EXNC + EX (LK) 

IF (LK.LE.NPLYO) TSTAR = TSTAR + T13(LK) 
IF (LK.GT.NPLYO) TSTAR = TSTAR + T12(LK) 


IF (LK.LE.NPLYO) ESTAR 
IF (LK.GT.NPLYO) ESTAR 


ESTAR + EX3(LK) 
ESTAR + EXX (LK) 


IF (l ZZ.EQ.2) TNMST = 0.0 
IF (I ZZ.EQ.2) GO TO 89 


TNMST = ( TNC - (TNC-TSTAR) *2*AL/WI DTH ) / ( EXNC - 
C (EXNC-ESTAR) *2*AL/WI DTH ) 

89 WRITE (6,*)' TNMST EQUALS ' .TNMST 


DO 164 LL = l.TPLY 

IF (LL.LE.NPLYO) WWC = (EX3(LL) * 

C (STRAIN(LST) + TNMST) - T13(LL) ) 

C * ( STRAIN(LST) - NMST3 (LL) + TNMST ) 

IF (LL.GT.NPLYO) WWC = (EXX(LL) * 

C (STRAIN(LST) + TNMST) -T12(LL) ) 

C * ( STRAIN(LST) - NMST2 (LL) + TNMST ) 

WWO = ( EX (LL) * (STRAIN(LST) + TNMST) - T1 (LL) ) 
C * ( STRAIN(LST) - NMSTO (LL) + TNMST ) 


164 GLC(JMM) = GLC(JMM) + WWO - WWC 


GLC(JMM) = GLC(JMM) / 2.0 

************************** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 



THIS IS TO CALCULATE THE INTERLAMINAR SHEAR STRESSES. 

************************************************************* ******** 
UNCL - (WIDTH / 2.0)- AL 
DO 180 JX - 80,100 
JY «( 1.0 - JX /100.0) * UNCL 

SIGX(JMM.JX) - 0 
SI GY (JMM, JX) « 0 
DO 180 JS = 1,4 

SIGX(JMM.JX) «' BG(JS) * S (JS) * DEXP ( -S (JS) * JY ) 

C * NXY 1 (JS) + SIGX (JMM, JX) 

S I GY (JMM, JX) = BG (JS) * S (JS) * DEXP ( -S (JS) * JY ) 

C * NYl(JS) + SI GY (JMM, JX) 

180 CONTINUE 


THE PROGRAM CONTINUES AND FINDS THE VARIOUS STRAIN ENERGY RELEASE 
COMPONENTS 


IF (COUPL.EQ.l) GO TO I 65 

THIS IS FOR A SYSTEM THAT IS COUPLED, THE CRACK LENGTH IS 

DEL - S(4) * S (2) * ATH * ATH 
DEL - DEL * DEL * 0.6144 
GO TO 170 

THIS IS FOR AN UNCOUPLED SYSTEM 

I 65 SSW = .65 * ( S(l) + S( 2 ) + S ( 3 ) ) 

DEL - 18.7 * S(4) * SSW * ATH * ATH 
DEL = DEL * DEL / 571-00 
170 DEL - 135-7 * DEL * ATH 

C IF (LIL.EQ.1) WRITE (6,21 1) DEL 

FY - ZR 
FX - ZR 


DO 175 JP=1 ,4 

CONY- NY1 (JP) *BG (JP) * ( DEXP ( -S (JP) * DEL ) - 1 ) 
CONXY- NXY 1 (JP) * BG (JP) * ( DEXP ( -S (JP) *DEL)- 1 ) 
CCONXY - CCONXY + CONXY/S (JP) 

175 CCONY - CCONY + CONY/S (JP) 

*************************************************************** 

*************************************************************** 

FY - BG1 + CCONY / DEL 
FX - BG2 + CCONXY / DEL 

GII(IZZ.JM) = FY / 2.0 * THETV *STRAIN(LST) 

Gill (IZZ.JM) = FX / 2.0 * THETU * STRAIN (LST) 

DIFFG- Gl I (IZZ, JM) - G I I I (I ZZ, JM) 

CON = 2 

IF (DIFFG. GT. REAL (GLC (JMM))) CON-1 

IF (DIFFG. GT. REAL (GLC (JMM))) DEL = DEL * -9 

IF (DIFFG. GT. REAL (GLC (JMM))) WRITE( 6 ,*)' IT EXPLODES 


Gl (IZZ.JM) = GLC (JMM) -Gl I (IZZ, JM) -Gill (IZZ.JM) 



200 


CONTINUE 


RESULTS ARE PRINTED FOR EACH RUN. 


300 CONTINUE 

WRITE (6,266) STRAIN(LST) 

WRITE (6,267) 

WRITE ( 6 , 269 ) GLC(O) , G I (2.1) , G I l (2.1). GMI(2.1), 

C G I (2,1) /GLC (0) 

DO 350 l-l.MMC +1 

WRITE (6,268) CMOIST(I) ,GLC (I) , G I (1 , I ) , G I I (1 , I ) , 

C Gl I I (1, I) , Gl (1 , I) /GLC (I) 

350 CONTINUE 

WRITE (6,287) 

DO 36 O NS - 80,90,2 

360 WRITE (6,285)NS/100., S I GX (0 , NS) , ( S I GX (KL , NS) , KL«1 , 22 , 4) 

DO 365 NS « 91,100 

365 WR I TE ( 6 , 285) NS/ 100 . , S I GX ( 0 , NS) , ( S I GX (KL , NS) , KL-1 , 22 , 4) 

WRITE (6,286) 

DO 370 NS - 80,90,2 

370 WRITE(6,285) NS/100. ,SIGY(0, NS) ,( SIGY(KL, NS) , KL-1, 22, 4) 

DO 375 NS « 91.100 

375 WRITE(6,285) NS/ 100 . ,SIGY(0, NS) ,( SIGY(KL, NS) , KL«1,22,4) 


400 CONTINUE 

201 FORMAT (//,' THE WIOTH OF THE LAMINATE IS \F 8 . 5 ) 

287 FORMAT (/////, ' THESE ARE THE IN-PLANE INTERLAMINAR SHEAR 
C 'STRESSES -- SIGMA XY ',/,' THEY ARE FOUND AT VARIOUS.', 

C '• MOISTURE CONTENTS •',//,' Y LOCAT I ON • , 7X, ' MECH 0NLY',8X, 

C 'H-0.0' , 10X, 'H-0.2' ,10X, 'H-0.4' , 10X, 

C 'H-0.6' , 10X, 'H-0.8' , 10X, 'H-l .0' ,//) 

285 FORMAT (3X,F7.2,4X,7F15*8) 

202 FORMATC THE NUMBER OF LAMINATES ABOVE AND BELOW THE CRACK IS' 

C . 13 . 5 x. 13 ) 

204 FORMAT (///,' THE PLYS ARE INPUTTED FROM BOTTOM TO TOP',/, 

. C 'BUT THE PLY CHARACTERISTICS FROM TOP TO BOTTOM ARE ') 

206 FORMAT (//, ' FOR PLY ',15,' THE SUBLAMINATE HAS THESE PROPERTIES') 

205 FORMAT (//,' WITH THIS LAYUP, THE PLYS ARE COUPLED ',//) 

210 FORMAT (//,' WITH THIS LAYUP, THE PLYS ARE DECOUPLED ',//) 

286 F ORMAT (//, ' THESE ARE THE OUT-OF-PLANE INTERLAMINAR SHEAR ', 

C 'STRESSES — SIGMA YZ ',/,' THEY ARE FOUND AT VARIOUS', 

C ' MOISTURE CONTENTS ',//,' Y LOCAT I ON ' ,8X, ' MECH 0NLY',8X, 

C * H=0 . 0 ' , 1 OX , ' H=0 . 2 ' , 1 OX , ' H=0 . 4 ' , 10X, 

C 'H-0.6' , 10X, 'H-0.8' ,10X, 'H-l .0' ,//) 

208 FORMATC El AND E2 ARE (MSI) ' ,F8.4, 10X.F8.4) 

289 FORMAT (//,' THE LAMINA PLY CHARACTERISTICS INITIALLY ARE ',/) 

288 FORMAT (/,' S1NM AND S2NM ARE EQUAL TO ' , F 1 4 . 10, 4X , F 1 4 . 10) 

211 FORMAT (//,' THE CRACK LENGTH STEP SIZE IS ' , F 1 2 . 8) 

266 FORMAT ('O' T , ' THE STRAIN IS EQUAL TO ' , F 1 2 . 7 , /, 

C ' THE VALUES OF GT, Gl, Gil, AND GUI ARE IN IN-LB/IN/IN 
C ') 

267 FORMAT (/,3X, '% CMOIST' ,8X, 'GGG (WHITNEY) ' ,6X, 'Gl ' ,9X, 

C ' Gl I ' ,8X, 'Gl I I ' ,6X, ' Gl/G (W-T) ' ,//) 

269 FORMAT (/, ' MECH. ONLY ' , 3X,F 12.9.4 ( 2 X,F 1 1 . 7 ) ,/) 

268 FORMAT(5X,F8.3,3X,F12.9.4(2X,F11.7) ) 

215 FORMAT (////, ' THETA V. IS ' , F 15 . 10 , ' THETA U IS ' , F 24 . 1 9) 

216 FORMAT (//, ' NY IS '.F23.ll,' NXY IS ' , F23 - 1 8) 



207 FORMAT (/,' THE THICKNESS AND THETA VALUES ARE ' , F9- 6 .5X. F8 . 3) 
209 FORMAT ( 1 THE POISSON RATIO (1.2) IS ' .F10.5) 

211* FORMAT ( ' G OF (1-2), AND (3~1) ARE -MSI \ 2(F9-A,2X)) 

217 F0RMAT('0' ,////, 8X, 'THE FOUR CHARACTERISTIC VALUES ASSOCIATED' 

C ,/,8X, ' WITH THE 8 DEGREE POLYNOMIAL FOR THE COUPLED CASE ARE') 

218 FORMAT (//,5X. ' THETA ' ,£>X, ' NY ' ,8X, ' NXY ' ,//) 

219 FORMAT (///,' THE S VALUES OF THE MEMBRANE ARE ' , F 1 5 -5. 3X, F 15 -5) 

220 FORMAT(/,LX,F9.l4,3X,F9.2,3X,F9.2) 

221 FORMAT (/, ' S OF ',12,' IS EQUAL TO \F20.10) 

223 FORMAT (/, ' THE UNCOUPLED SY (1,2) VALUES ARE ' , F 15 -5 . 3X, F 1 5 -5) 

22A FORMAT (/,' THE UNCOUPLED SX (1,2) VALUES ARE * ,F 15-5.3X.F 15-5) 

231 FORMAT (//,' THE STRAIN IS EQUAL TO '.F12.8,/, 

C ' THE CHANGE IN TEMPERATURE IS '.F12.5,/, 

C ' THE COEFFICIENTS SWELLING DUE TO MOISTURE ARE ' , 2 (2X, F 1 2 .8) 

C THE COEFFICIENTS OF THERMAL EXPANSION ARE ' ,2 (2X,F I 5 . 9 ) ) 


232 FORMAT (/,' THE MOISTURE COEFF I CENT IS ' . F 1 5 . 8 ) 

233 FORMAT (/,' THE MOISTURE COEFFICIENT VARIES FROM 0 TO 1.2 ') 



THE LAMINA PLY CHARACTERISTICS INITIALLY ARE 
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THE WIDTH OF THE LAMINATE IS 1.51200 
THE NUMBER OF LAMINATES ABOVE AND BELOW THE CRACK IS 3 1 


FOR PLY 1 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .0051*00 35.000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1.2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 2 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE . 005^00 -35-000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 3 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .0051*00 
El AND E 2 ARE (MSI) I 8 . 7 OOO 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1)' ARE -MSI .8320 .8320 


.000 

1.2300 


FOR PLY A THE SUBLAMINATE HAS THESE PROPERT 
THE THICKNESS AND THETA VALUES ARE .0051*00 
El AND E2 ARE (MSI) 18.7000 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 


IES 

90.000 

1.2300 


.8320 


THE STRAIN IS EQUAL TO .00251*000 
THE CHANGE IN TEMPERATURE IS -280.00000 

THE COEFFICIENTS SWELLING DUE TO MOISTURE ARE .00000000 

THE COEFFICIENTS OF THERMAL EXPANSION ARE -.000000230 

THE MOISTURE COEFFICIENT VARIES FROM 0 TO 1.2 


WITH THIS LAYUP, THE PLYS ARE COUPLED 


S 

S 

S 

S 


THE FOUR CHARACTERISTIC- VALUES ASSOCIATED 
WITH THE 8 DEGREE POLYNOMIAL FOR THE COUPLED CASE ARE 


OF 

1 

IS 

EQUAL TO 

1 * 07 . 057368271 * 1 * 

OF 

2 

IS 

EQUAL TO 

11 * 1.1 1977801*18 

OF 

3 

IS 

EQUAL TO 

1 16.7332723860 

OF 

1 * 

IS 

EQUAL TO 

55.551*1*207729 


.00556000 

. 000011*900 


THE UNCOUPLED SX (1,2) VALUES ARE 


392.221*78 


106.21737 


THE UNCOUPLED SY (1,2) VALUES ARE 


13^.68589 


56 . 2575 ** 


THE S VALUES OF THE MEMBRANE ARE 1 77 - 389^5 66.26466 


S1NM AND S2NM ARE EQUAL TO -.0000157292 .0000000000 

SI AND S2 = 1 34.932471163957803943120257 641.500299099584182787943089 
SX AND SY = 7.OO3582O8142O9O671392409983 33.2966554398959572557143052 


THETA V IS .3932559441 THETA U IS -.0981017742335529623 
NY IS -38.32041690358 NXY IS -61.959633961712929972 


THE STRAIN IS 

EQUAL TO .0025400 

THE VALUES OF 

GT, Gl, Gil, 

AND GUI 

% CMOIST 

GGG (WHITNEY) 

Gl 

MECH. ONLY 

.101653935 

.0670191 

.000 

•522939955 

.4084740 

.050 

.488576635 

•3794622 

.100 

.455156902 

.3513864 

.150 

.422680755 

.3242463 

.200 

. 39"S 148194 

.2980421 

.250 

.360559220 

.2727737 

.300 

.330913832 

.2484412 

• 350 

.302212031 

.2250445 

.400 

.274453816 

.2025836 

.450 

.247639187 

.1810586 

.500 

.221768145 

.1604694 

.550 

.196840689 

.1408161 

.600 

.172856820 

.1220986 

.650 

.149816537 

.1043169 

• 700 

.127719840 

.0874710 

.750 

. 106566730 

.0715610 

.800 

.086357206 

.0565869 

.850 

.067091269 

.0425486 

.900 

.048768918 

.0294461 

• 950 

.031390154 

.0172794 

1 .000 

.014954976 

.0060486 

1.050 

-.000536616 

-.0042464 

1.100 

-.015084621 

-.0136055 

1.150 

-.028689040 

-.0220288 

1 .200 

-.041349872 

-.0295163 


ARE IN IN-LB/IN/IN 


Gl 1 

Gl 1 1 

Gl/G (W-T) 

.0346174 

.0000174 

.6592868 

.1144151 

.0000509 

.7811106 

. 1090658 

.0000486 

•7766688 

.1037241 

.0000464 

.7720115 

.0983903 

.0000442 

.7671187 

.0930641 

.0000420 

.7619672 

.0877458 

.0000397 

.7565296 

.0824351 

.0000375 

.7507731 

.0771323 

.0000353 

.7446576 

.0718371 

.0000331 

.7381338 

.0665497 

.0000308 

.7311387 

.0612701 

.0000286 

.7235909 

.0559982 

.0000264 

.7153809 

.0507341 

.0000242 

.7063566 

»0454777 

.0000220 

.6962975 

.0402291 

.0000197 

.6848665 

.0349882 

.0000175 

.6715139 

.0297550 

.0000153 

.6552653 

.0245296 

.0000131 

.6341891 

.0193120 

.0000109 

.6037875 

.0141021 

.0000087 

.5504724 

.0088999 

.0000064 

.4044538 

.0037055 

.0000042 

7.9132595 

-.0014811 

.0000020 

.9019465 

-.0066600 

-.0000002 

.7678481 

-.0118312 

-.0000024 

.7138182 



THESE ARE THE IN-PLANE INTERLAMINAR SHEAR STRESSES -- SIGMA XZ 
•l HEY ARE FOUND AT VARIOUS MOISTURE CONTENTS 
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THE LAMINA PLY CHARACTERISTICS INITIALLY ARE 


THE WIDTH OF THE LAMINATE IS 1.51200 
THE NUMBER OF LAMINATES ABOVE AND BELOW THE CRACK IS 3 1 

THE PLYS ARE INPUTTED FROM BOTTOM TO TOP 
BUT THE PLY CHARACTERISTICS FROM TOP TO BOTTOM ARE 


FOR PLY 1 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .0051*00 35.000 

El AND E2 ARE (MSI) 18-7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2). AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 2 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE . 005^00 

El AND E2 ARE (MSI) 18-7000 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


.000 

1.2300 


FOR PLY 3 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005^00 -35. 

El AND E2 ARE (MSI) 18.7000 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


000 

1.2300 


FOR PLY It THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005^00 90. 

El AND E2 ARE (MSI) • 18.7000 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


000 

1.2300 


THE STRAIN IS EQUAL TO .00251*000 
THE CHANGE IN TEMPERATURE IS -280.00000 

THE COEFFICIENTS SWELLING DUE TO MOISTURE ARE .00000000 

THE COEFFICIENTS OF THERMAL EXPANSION ARE -.000000230 

THE MOISTURE COEFFICIENT VARIES FROM 0 TO 1.2 


WITH THIS LAYUP, THE PLYS ARE COUPLED 


THE FOUR CHARACTERISTIC VALUES ASSOCIATED 
WITH THE 8 DEGREE POLYNOMIAL FOR THE COUPLED CASE ARE 


OF 

1 

IS 

EQUAL 

TO 

360. 71 621*235^3 

OF 

2 

IS 

EQUAL 

TO 

136.396160LL92 

OF 

3 

IS 

EQUAL 

TO 

113-9856581*772 

OF 

*4 

IS 

EQUAL 

TO' 

5b. 080 1738692 


.00556000 

. 0000 U 900 



THE UNCOUPLED 

SX (1,2) VALUES 

ARE 

352.6287*4 

86.89277 


THE UNCOUPLED 

SY (1,2) VALUES 

ARE 

126. *45615 

59.53405 

63 

THE S VALUES OF 

THE MEMBRANE ARE 193-07807 

70.65819 


S1NM AND S2NM 

ARE EQUAL TO 

■.0000157292 .0000000000 


SI AND S2 

= 1314.932*4711639578039^3120257 614 1 .50029909958*4 1827879*43089 

SX AND SY 

■= 7.003582081142090671392*409983 33.296655*43989595725571*43052 

THETA V IS 

• 978*402*4562 THETA U IS 

.00000000000000097 1 8 


NY IS 

-38.320*4 1690358 

NXY IS - 

■61.959633961712929972 


THE STRAIN IS 

EQUAL TO .0025*400 




THE VALUES OF 

GT, Gl, Gil, 

AND Gill 

ARE IN IN- 

LB/IN/IN 


% CMO 1 ST 

GGG (WHITNEY) 

Gl 

Gl 1 

Gill 

Gl/G (W-T) 

MECH. ONLY 

.09*4207177 

.0078275 

.0863797 

.0000000 

.0830882 

.000 

.510289070 

.2262277 

.284061*4 

.0000000 

.4433324 

.050 

.*476*405835 

.2055316 

.2708742 

.0000000 

.4314212 

. TOO 

.*♦143*4*46720 

. 1857*497 

.2576970 

.0000000 

.4188771 

.150 

.*411*411727 

. 1668819 

.2445298 

.0000000 

.4056324 

.200 

.380300856 

. 1*48928*4 

.2313725 

.0000000 

.3916067 

.250 

.3501 1*4106 

.1318890 

.2182251 

.0000000 

.3767028 

.300 

.320851*477 

•1157637 

.2050878 

.0000000 

.3608016 

• 350 

.292512969 

.1005526 

.1919603 

.0000000 

.3437545 

. 400 

.265098582 

.0862557 

.1788428 

.0000000 

.3253723 

.1*50 

.238608317 

.0728730 

•1657353 

.0000000 

.3054085 

.500 

.2130*42173 

.060*40*4*4 

•1526377 

.0000000 

.2835328 

.550 

. 188*400151 

.01488500 

.1395501 

.0000000 

.2592888 

.600 

. 16*4682250 

.0382098 

.1264724 

.0000000 

.2320215 

.650 

. 1*41888*470 

.028*4838 

.1134047 

.0000000 

.2007475 

.700 

.12001881 l 

.0196719 

.1003469 

.0000000 

.1639065 

.750 

.09907327*4 

.01177*41 

.0872991 

.0000000 

.1188427 

.800 

.079051858 

.00*47906 

.0742613 

.0000000 

.0606005 

.850 

.05995*4563 

-.0012788 

.0612334 

.0000000 

-.0213295 

.900 

.0*4 1 781389 

-.006*43*40 

.0482154 

.0000000 

-.1539925 

.950 

.02*4532337 

-.0106751 

.0352074 

.0000000 

-.435142° 

1 .000 

.008207*406 

-.01*40020 

.0222094 

.0000000 

-1.7060144 

1.050 

-.007193*403 

-.016*41*47 

.0092213 

.0000000 

2.2819054 

1 . 100 

-.021670092 

-.0179132 

-.0037569 

.0000000 

.8266329 

1.150 

-.035222659 

-.018*4976 

-.0167251 

.0000000 

.5251615 

1 .200 

-.0*47851 10*4 

-.0181678 

-.0296833 

.0000000 

.3796734 


THESE ARE THE IN-PLANE INTERLAMINAR SHEAR STRESSES -- SIGMA XZ 
THEY ARE FOUND AT VARIOUS MOISTURE CONTENTS 
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THE LAMINA PLY CHARACTERISTICS INITIALLY ARE 
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THE WIDTH OF THE LAMINATE IS 1.51200 
THE NUMBER OF LAMINATES ABOVE AND BELOW THE CRACK IS 2 


FOR PLY 1 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 30.000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 2 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .0051*00 -60.000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 3 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 75-000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 4 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 -15-000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


THE STRAIN IS EQUAL TO .00251*000 
THE CHANGE IN TEMPERATURE IS -280.00000 

THE COEFFICIENTS SWELLING DUE TO MOISTURE ARE .00000000 .00556000 

THE COEFFICIENTS OF THERMAL EXPANSION ARE -.000000230 .000014900 

THE MOISTURE COEFFICIENT VARIES FROM 0 TO 1.2 


WITH THIS LAYUP, THE PLYS ARE COUPLED 


THE FOUR CHARACTERISTIC VALUES ASSOCIATED 
WITH THE 8 DEGREE -POLYNOMIAL FOR THE COUPLED CASE ARE 

S OF IIS EQUAL TO 202 . 3666962066 

S OF 2 IS EQUAL TO 150.2447231*209 

S OF 3 IS EQUAL TO 96.3990023366 

S OF 4 IS EQUAL TO 72.1855545293 


THE UNCOUPLED S>: (1.2) VALUES ARE 


178.08581 


9C.73315 



THE UNCOUPLED SY (1,2) VALUES ARE 


137.79321. 


77. 57481 


THE S VALUES OF THE MEMBRANE ARE 107-30161. 75-55664 

S1NM AND S2NM ARE EQUAL TO -.0013126127 .0012407192 

SI AND S2 « li«5 -7 12738391*363 11*5201 23A 301* 261.71259626539561223751*2175 
SX AND SY ■= 35.5756851 1671*391 1519011*7879 60.958551*9677311*312319553312 


THETA V IS .2402517909 THETA U IS 2 . 07 31*381+470 2 1 801786 1 
NY IS -35.1*8501616832 NXY IS -6l.l*6l8509109l*0726086 


THE STRAIN IS EQUAL TO .00251*00 


THE VALUES OF 

GT, Gl, Gil, 

AND Gill 

ARE IN IN- 

LB/IN/IN 


S CMOIST 

GGG (WHITNEY) 

Gl 

G 1 1 

Gl 1 1 

G 1 /G (W-T) 

MECH. ONLY 

.171*065036 

.0136731* 

.0093651 

.1510265 

.0785535 

.000 

.28859921*1* 

.0062276 

.0365562 

• 21*58155 

.0215787 

.050 

.277169512 

.0030909 

.031*71*33 

.2393353 

.011 1516 

.100 

.266289511* 

.0005010 

.0329323 

.2328563 

.0018813 

.150 

.255959252 

-.00151*22 

.0311230 

.2263784 

-.0060250 

.200 

.21.6178721* 

-.0030385 

.0293158 

.2199016 

-.0123427 

.250 

.23691*7932 

-.0039880 

.0275100 

.2131*260 

- .0168309 

.300 

.228266871* 

-.001*3908 

.0257062 

.2069515 

-.0192354 

• 350 

•220135552 

-.001*21*68 

.023901*2 

.2004781 

-.0192916 

.1*00 

.212553961* 

-.0035559 

.022101*0 

.1940059 ‘ 

-.0167296 

.1*50 

.205522112 

-.0023183 

.0203057 

.1875347 

-.0112801 

.500 

.199039991* 

-.0005339 

.0185092 

.1810647 

-.0026824 

• 550 

.193107612 

.0017973 

.016711*1* 

.1745959 

.0093072 

.600 

. 187721*965 

.001*6753 

.011*9215 

.1681281 

.0249050 

.650 

. 182892052 

.0081001 

.0131301* 

.1616615 

.0442888 

.700 

. 178608875 

.0120716 

.01131*12 

.1551961 

.0675870 

• 750 

.17i*875i*32 

.0165900 

.0095537 

.1487317 

.0948676 

.800 

.171691725 

.0216552 

.0077681 

. 1422685 

.1261282 

.850 

.169057753 

.0272671 

.005981*2 

.1358064 

. 1612888 

• 900 

.166973515 

.0331*259 

.001*2022 

.12931*51* . 

.2001866 

• 950 

.1651*39013 

.01*01311* 

.0021*220 

.1228856 

.2425751 

1 .000 

.161*1*51*21*5 

.01*73837 

.0006L36 

.1164269 

.2881270 

1 .050 

.161*019213 

.0551828 

-.0011329 

.1099693 

.3364412 

1 . 100 

. 161*133916 

.0635287 

-.0029077 

.1035129 

.3870542 

1.150 

.161*798353 

.0721*211* 

-.001*6806 

.0970576 

.4394548 

1.200 

.166012526 

.0818609 

-.0061*517 

.0906034 

.4931008 



THESE ARE THE IN-PLANE INTERLAMINAR SHEAR STRESSES -- SIGMA XZ 
THEY ARE FOUND AT VARIOUS MOISTURE CONTENTS 
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THE LAMINA PLY CHARACTERISTICS INITIALLY ARE 


THE WIDTH OF THE LAMINATE IS 1-51200 88 

THE NUMBER OF LAMINATES ABOVE AND BELOW THE CRACK IS 2 2 


FOR PLY 1 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .0051*00 -35-000 

El AND E2 ARE (MSI) 18-7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 2 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 55-000 

El AND E2 ARE (MSI) 1 8 . 7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 3 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 10.000 

El AND E2 ARE (MSI) 1 8. 7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


FOR PLY 4 THE SUBLAMINATE HAS THESE PROPERTIES 
THE THICKNESS AND THETA VALUES ARE .005400 -80.000 

El AND E2 ARE (MSI) 18.7000 1.2300 

THE POISSON RATIO (1,2) IS .29200 

G OF (1-2), AND (3-1) ARE -MSI .8320 .8320 


THE STRAIN IS EQUAL TO .00254000 
THE CHANGE IN TEMPERATURE IS -280. 00000 

THE COEFFICIENTS SWELLING DUE TO MOISTURE ARE 
THE COEFFICIENTS OF THERMAL EXPANSION ARE 


.00000000 

-.000000230 


.00556000 

.000014900 


THE MOISTURE COEFFICIENT VARIES FROM 0 TO 1.2 


WITH THIS LAYUP, THE PLYS ARE COUPLED 


THE FOUR CHARACTERISTIC VALUES ASSOCIATED 
WITH THE 8 DEGREE POLYNOMIAL FOR THE COUPLED CASE ARE 


S OF 

1 

IS EQUAL 

TO 

233.9388572236 


S OF 

2 

IS EQUAL 

TO 

156.9014619788 


S OF 

3 

IS EQUAL 

TO 

115.2992565645 


S OF 

4 

IS EQUAL 

TO 

48.0575100718 


THE 

UNCOUPLED SX 

(1.2) 

VALUES ARE 186.41344 

96.42124 

THE 

uncouple:- S'- 


VALUES ARE 128.78572 

49.37034 



THE S VALUES OF THE MEMBRANE ARE 


1 19.36060 


57.91987 


89 


S1NM ANO S2NM ARE EQUAL TO - .0007221 149 -.0007139234 

SI AND S2 - H3.69l863395728751259l*52089 285.261.127101.852281 15901*5526 
SX AND SY - 32. 1436979546791447820272102 62.201.9978939766992050711339 


THETA V IS .1*31.7536621 THETA U IS -1 .731*44 151.08201 789576 
NY IS -51.. 366 1988991.6 NXY IS 1.5.6186571.1.5050692111 


THE STRAIN IS EQUAL TO .00251.00 


THE VALUES OF 

GT, Gl, Gil, 

AND Gill 

ARE 4 N IN- 

•LB/IN/IN 


% CMOIST 

GGG (WHITNEY) 

Gl 

G 1 1 

Gill 

G 1 /G (W-T) 

MECH. ONLY 

.131525889 

.0207693 

.0219689 

.0887877 

.1579102 

.000 

.391.637005 

.1683076 

.0866210 

.1397081* 

.1*261*871 

.050 

.36931.7077 

.1508205 

.0823111* 

.1362152 

.40831*35 

.100 

.3^5181*223 

.131*1*556 

.0780061 

.1327225 

.3895185 

. 150 

.32211*81*1*3 

.1192131 

.0737050 

. 1292301* 

.3700561* 

.200 

.300239735 

. 1050929 

.0691*082 

.1257387 

.3500298 

.250 

.2791*58102 

.092091*9 

.0651156 

. 12221*75 

.3295483 

.300 

.2598035lt2 

.0802193 

.0608271* 

.1187569 

.3087691 

•350 

.21*1276055 

.0691.660 

.05651*31* 

.1152667 

.2879108 

.1.00 

.22387561*2 

.0598350 

.0522636 

•1117771 

.2672688 

.1*50 

.207602302 

.0513262 

.01*79882 

. 1082879 

.2472335 

.500 

.1921*56036 

.01*39398 

.01*37170 

. 101*7992 

.22831 10 

.550 

. 1781*3681*3 

.0376757 

.0391*500 

.1013111 

.21 1 1431 

.600 

. l655l*l*721* 

.0325339 

.0351871* 

.0978235 

.1965261* 

.650 

.153779678 

.028511*1* 

.0309290 

.091.3363 

.1854237 

.700 

. 1 1*3 1 1+ 1 706 

.0256172 

.026671*8 

.09081*97 

.1789638 

.750 

.133630807 

.02381*23 

.0221*250 

.0873635 

.1784190 

.800 

.12521*6982 

.0231897 

.0181791* 

.0838779 

.1851516 

.850 

.117990230 

.0236591* 

.0139381 

.0803928 

.2005198 

.900 

. 1 11860552 

.0252511* 

.0097010 

.0769081 

.2257398 

•950 

. 10685791*7 

.0279657 

.0051*683 

.073421.0 

.2617089 

1 .000 

. 1029821*16 

.0318023 

.0012397 

.06991*01* 

.3088128 

1 .050 

.100233958 

.0367612 

-.002981*5 

.0661*573 

•3667539 

1.100 

.098612573 

• Ol.281.2l. 

-.007201*5 

.0629747 

.4344517 

1.150 

.0981 18262 

.05001*59 

-.01 14202 

.0594925 

.5100571 

1 .200 

.098751025 

.0583717 

-.0156316 

.0560109 

.5910999 



THESE ARE THE IN-PLANE INTERLAMINAR SHEAR STRESSES -- SIGMA X2 
THEY ARE FOUND AT VARIOUS MOISTURE CONTENTS 
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FRACTURE ANALYSIS OF LOCAL OELAMINATIONS IN LAMINATED COMPOSITES 

P. Sri ram and E. A. Armani os 
School of Aerospace Engineering 
Georgia Institute of Technology 
Atlanta, Georgia 30332 

Abstract 

Delamination is a predominant failure mode in continuous fiber rein- 
forced laminated composite structures. One type of delamination is the 
transverse crack tip delamination which originates at the tip of transverse 
matrix cracks. An analytical model based on the sublaminate approach and 
fracture mechanics is developed in this paper to study the growth of such 
delaminations. Plane strain conditions are assumed and estimates are 
provided for the total strain energy release rate as well as the mode 1 and 
mode II contributions. The energy release rate estimates are used to 
predict critical delamination growth strains and stresses by assuming a 
critical energy release rate. These predictions are compared with experi- 
mental data on T30Q/934 Graphite Epoxy [±25/90 n ] s laminates in the range 
n=.5 to 8. A good agreement is demonstrated for the range of n where the 
experimental observations indicate transverse crack tip delamination to be 
the predominant failure mode. 

Introduction 

Fiber reinforced composites are now being used in a wide variety of 
engineering structures. The concept of directional strength and stiffness 
has been, for the most part, understood sufficiently to enable efficient 
load bearing designs. One of the current major issues in composite struc- 
tures is the understanding and prediction of damage modes and failure 
mechanisms. A thorough knowledge of the failure mechanisms is oound to 
lead to the design of efficient and durable structures. Failures in these 
materials often initiate in the form of matrix cracks or delaminations. 


Matrix cracks refer to intralaminar failures whereas delaminations refer to 


interlaminar failures. 

Matrix cracks usually occur within laminates where the fibers run at 
an angle to the primary load direction. Hence, such matrix cracks are also 
called transverse cracks. Based on the location and direction of growth, 
two distinct types of delamination can be discerned. These two types are 
called edge delamination and local or transverse crack tip delamination. 
Edge delaminations initiate at the load free edges of the structure whereas 
local delaminations start from a transverse matrix crack. In many cases, 
both types occur concurrently with varying levels of interaction. It has 
been observed in simple tension tests of uniform rectangular cross section 
specimen (Edge Delamination test) that delaminations initiate along the 
load free edges and propagate normal to the load direction. Transverse 
matrix cracks running parallel to the fibers have also been observed in off 
axis plies such as 90° plies. Such transverse cracks terminate where the 
ply orientation changes. Delaminations can originate at the interface 
where transverse cracks terminate. These delaminations, called transverse 
crack delaminations or local delaminations, grow normal to the transverse 
crack from which they originate. In the case of 90° plies, the growth 
direction is parallel to the load. 

The growth process of edge delaminations and local delaminations is 
often modelled using a fracture mechanics approach leading to the calcula- 
tion of a strain energy release rate. This is because the strain energy 
release rate can correlate delamination behavior from different loading 
conditions and can account for geometric dependencies. The strain energy 
release rate associated with a particular growth configuration is a measure 
of the driving force behind that failure mode. in combination with 
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appropriate failure criteria, the strain energy release rate provides a 

means of predicting the failure loads of the structure. 

Several methods are available in the literature for analyzing edge 

1-3 

del ami nations. These include finite element modelling , complex variable 

4 

stress potential approach , simple classical laminate theory based tech- 
nique and higher order laminate theory including shear deformations . 
Finite element models provide accurate solutions but involve intensive 
computational effort. Classical laminate theory (CLT) based techniques 
provide simple closed form solutions and are thus well suited for prelimi- 
nary design evaluation. Classical laminate theory based techniques provide 
only the total energy release rate, and thus in a mixed mode situation, 
there is insufficient information to completely assess the delamination 
growth tendency. A higher order laminate theory including shear deforma- 
tions has the ability to provide the individual contributions of the three 
fracture modes while retaining the simplicity of a closed form solution. A 
shear deformation model is available for edge delamination and has been 

/r 

shown to agree well with finite element predictions . 

Crossman and Wang^ have tested T300/934 Graphite epoxy [±25/90^3 s 
specimens in simple tension and reported a range of behavior including 

O 

transverse cracking, edge delamination and local delamination. O'Brien 
has presented classical laminate theory solutions for these specimen, 
demonstrating reasonabl e agreement in the case of edge del ami nation but 
with some discrepancies in the local delamination predictions. An empiri- 
cal finite element based combined edge and local delamination formulation 

g 

has also been proposed . Its predictions, however, do not fully explain 
the dependency of the critical strain on the number of 90° plies. 

In this paper, a shear deformation model is developed for the analysis 
of local delaminations originating from transverse cracks in 90° plies 


located in and around the specimen midplane. Plane strain conditions are 
assumed and thickness strain is neglected. Delaminations are assumed to 
grow from both ends of the transverse crack tip. The transverse crack is 
treated as a free boundary and the delamination is considered to be the 
crack whose growth behavior is to be modelled. The sublaminate ap- 
proach 1 ^’ 11 is used to model different regions of the specimen. The 

resulting boundary value problem is solved to obtain the interlaminar 
stresses, total strain energy release rate and energy release rate compo- 
nents. Critical local delamination growth loads are predicted for the 
[±25/90 n ] s specimen. 

Analytical Model 

The formulation is based on the sublaminate approach detailed in ref. 
10. A longitudinal section illustrating the geometry of a generic configu- 
ration is shown in fig. 1. The central region is assumed to be made of 90° 
plies with an isolated transverse crack in the middle. Oelaminations are 
assumed to grow from both ends of the transverse crack, and towards both 
ends as shown. From symmetry considerations, only one quarter of the 
configuration is modelled. The modelled portion is divided into four 
sublaminates as shown in fig. 2. The top surface (subiaminates 1 and 4) is 
stress free. In order to simplify the analysis, plane strain conditions 
are assumed and the thickness strain (e z ) is set to zero. The consequence 
of this combined with the fact that the w displacement is zero along the 
center line is that w is zero in sublaminates 1,2 and 3. Further, this 
approximation does not allow for the enforcement of boundary conditions on 
the shear stress resultants, leading to incorrect estimates of the inter- 
laminar normal stresses. The interlaminar shear stresses, however, are not 

fi in 

affected by this assumption ’ . The assumptions lead to considerable 

simplifications in the analysis. In spite of the simplifications, reliable 


4 


energy release rate components can be estimated based on the interlaminar 
shear stress distributions^’*^ 1 . 

A generic sublaminate is shown in fig. 3 along with the notations and 
sign conventions. The peel and interlaminar shear stresses are denoted by 
P and T respectively with t and b subscripts for the top and bottom surface 
respectively. The axial stress resultant, shear stress resultant and 
bending moment resultant are denoted by N, Q and M respectively. A summary 
of the governing equations is presented here for convenience. These are 
derived for a generic sublaminate using the principle of virtual work in 
Reference 12. 

The x and z displacements within the sublaminate are assumed to be of 
the form 

u(x,z)=U(x)+zp(z) (1) 

w(x,z)=W(x) . (2) 

Here U represents the axial midplane stretching and W is the transverse 
displacement. The shear deformation is recognized through the rotation 6. 
The origin of the coordinate axes for the sublaminates is taken at the 
delamination tip as shown in fig. 4. The equilibrium equations take the 
form 


N ,x + VV° 

(3) 

W p b=° 

(4) 

M x -Q+(h/2)(T t +T b )=0. 

(5) 


where h is the thickness of the sublaminate. The constitutive relations in 


terms of the force and moment resultants are 



where the A.., B. . and B-. are the classical laminate theory axial, cou- 
^ J ^ J ^ J 

pling and bending stiffnesses. The boundary variables to be prescribed at 
the sublaminate edges are 

N or U 
M or p 
Q or W. 

Additionally, at the interfaces between sublaminates, reciprocal traction 
and displacement matching boundary conditions have to specified. 

Solution Procedure 

A detailed solution is provided in the Appendix. A brief summary is 

provided here for convenience. The variables in sublaminates i and 2 are 

coupled by their reciprocal interlaminar stresses denoted T^ and P 1 and by 

displacement continuity at their common interface. Assuming exponential 

sx 

solutions for the axial force and bending moment resultants (N^=Ae , 
M 1 =Be sx etc.) leads to an eigen value problem involving the parameter s. 
The eigen values turn out to be 0 and two nonzero values (say s 1 and s 2 ) 
occurring in positive and negative pairs. Since the resultants maintain 
finite values as x tends to large negative values (left end of sublaminates 
1 and 2), the negative roots are dropped out of the solution. 

The following boundary conditions from the ends of the modelled region 
are enforced. 


n 2 (0)=0 

(9) 

Q 4 (a)=0 

(10) 

P 4 (a)=0 

(ID 

N l +N 2 = APPl i e d L°ad 

(12) 


Further, the following displacement matching conditions are applied. 

u 1 (x,-.5h 1 )=u 2 (x,.5h 2 ) 

u 1 (0)=u 4 (0) 
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(13) 

(14) 


u 2 ( 0 )=u 3 ( 0 ) 

(15) 

P 1 (0)=p 4 (0) 

(16) 


It should be noted that a ^ an< ( P 3 matching condition cannot be applied at 

this level of modeling since it would amount to specifying both W and 
6 12 

Q ’ . Consequently, there is a displacement discontinuity at the de- 
lamination tip. The effect of this will be discussed subsequently. To 

eliminate rigid body displacements, is set to zero at the left end. The 
following solutions can then be obtained for the resultants in sublaminates 
1 and 2 . 

N 1 =a 1 e s l x +a 2 e s 2 x +eA n(1) (17) 

N 2 =-a 1 e s l x -a 2 e s 2 x +eA n(2) (18) 

Ml= a i k i eSlX + a 2 k 2e S 2 x (19) 

M 2 =a 1 k 3 e s l x +a 2 k 4 e s 2 x ( 20 ) 

The interlaminar shear and peel stresses between sublaminates 1 and 2 can 
be obtained as 

T l= a l s ie SlX + a 2 s 2e S 2 X (21) 

P 1 =(k 1 +.5h 1 )(a 1 s 1 2 e s l x )+(k 2 +.5h 1 )(a 2 s 2 2 e s 2 x ) ( 22 ) 

In the above solutions, the k parameters are dependent on the eigen values 
and the stiffness of sublaminates 1 and 2 , the a parameters depend on the k 
parameters and the initial crack length a, and e is defined as 

c-o ( h ^+h 2 ) / ( A^^ ^ j+A^^ 2 ^ ) (23) 

where o is the applied uniform axial stress. Complete expressions for the 

eigen values and the a and k parameters can be found in the Appendix. 

Proceeding on to sublaminates 3 and 4, the following solutions can be 
written. 



V° 

(24) 


M 3 =<t>^sinh ^x+^^cosh w 3 x 

(25) 

where 

< * l 2'" a l* < '3 +a 2* < 4 ’ 

(26) 
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and 


^^“** 1^2 coth uj^ a 


0.5 


w 3 =(A 55(2) /D 11(2) ) 
N 4 =e(A ll(l) +A ll(2) ) 
M 4 =a l k l +a 2 k 2 


(27) 

(28) 

(29) 

(30) 


The corresponding displacement solutions are provided in the Appendix. 
The compliance of the specimen can be evaluated as 

C=2U 4 (a)/P (31) 

where P/2 is the load applied to the modelled section. The total energy 
release rate for the modelled section i.e. the total energy release rate 
per crack is then given by 

G t =P 2 /2 w (dC/da) (32) 

where w is the specimen width. Use of the previously described solutions 


leads to the following expression. 




(33) 


2w J (*4u(i) -Aii(i) + ^ii(2) + I' 2 

where the quantities and contain exponential terms dependent on the 
initial delamination length. Using the virtual crack closure technique, 
from the relative displacements in the cracked portion and the interlaminar 
stresses ahead of the crack tip, the mode I and mode II energy release rate 
contributions can be obtained. The mode III energy release rate is zero 
from the assumption of plane strain. The mode II energy release rate is 
given by 

n i: L [ t< ~ 

(34) 


G// = Um^ J T\(x — 6)t±u(x)dx 


where 6 is the virtual crack step size. The result of the limiting process 
is zero if there is no singularity in the stress field 10 . So, the limit is 
usually taken as the crack step size 6 tends to a small value, say a, based 
on the decay length or the length required to capture the essential fea- 
tures of the stress and displacement fields near the crack tip. The decay 
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In this study, the 


length is dependent on the eigen values s^ and s^. 
value of A has been set to 

Aq=. 25( 1/Sj + l/s 2 ) (35) 
since it reasonably fulfills the criterion given above. In a similar 
fashion, the mode I energy release rate can be obtained based on the normal 
stress (P) and the w displacements near the crack front. The normal (peel) 
stress estimate is inaccurate due to the absence of thickness strain. 
Hence, an alternate approach was used to estimate Gj, the mode I energy 
release rate. The total energy release rate for this problem is made up 
entirely of Gj and Gjj (Gjjj= 0). From an estimate of Gy and Gjj, an 
estimate for Gj can be obtained simply as 

G r G T ‘G II (36) 
The critical load for a given specimen, can then be evaluated based on an 
appropriate fracture law. This is illustrated in the following section. 

Results and Discussion 

The solutions derived in the previous section have been used to model 
the behavior of [±25/90 n ] s T300/934 Graphite Epoxy specimen for n values of 
.5, 1,2, 3, 4, 6, and 8. These correspond to the specimen tested by Crossman 
and Wang^. The specimen width and length were fixed at .0381 m and .015m 
respectively, as in the tests. The solutions were generated using a simple 
computer program based on the closed form expressions for the interlaminar 
stress and energy release rates. The applied load was set to 100 MPa, of 
the same order as in the tests. 

An example of the total energy release rate variation with the cracK 
length is presented in fig. 5. The asymptotic value of Gy is denoted by 
Gy 0 in the figure. It can be observed that after a certain crack length, 
the Gy is independent of the crack length. On the basis of curves like the 
one shown in fig. 5, the crack length was fixed at 10 ply thicknesses for 
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the remainder of the study. The dependence of the mode II contribution of 
the energy release rate on initial crack length (a) is depicted in fig. 6. 
Typical interlaminar shear and normal stress profiles are presented in 
figs. 7 and 8 respectively. The corresponding energy release rates have 
also been calculated and are presented in Table 1 and fig. 9. 

In order to evaluate the critical loads, an appropriate mixed mode 
fracture law has to be applied, based on the calculated energy release 
components. Since the calculated mode split shows only a small variation 
with n, the simple Griffith criterion G-j-=G-]- c has been used to scale the 

stresses to obtain the critical delamination growth stress (o c ) and strain 

2 

(e c ) values. The critical energy release rate Gy c was chosen as 415 J/m 
to obtain the critical stresses and strains listed in Table 1. This value 
of G^ c is larger than Gj c to account for the presence of mode II and the 
fact that G tt is about four times G T for the material system under 
consideration. The critical strains are plotted against n, the number of 
90° plies in fig. 10. The experimental results of ref. 7 and the predic- 
tions of refs. 8 and 9 are also presented in the figure for comparison. 
The predictions of the model developed in this paper are represented by the 
solid line while the experimental results are shown as filled squares. The 
classical laminate theory and finite element critical strain predictions of 
refs. 8 and 9 are represented by triangles with a connecting line and a 
dotted line respectively. 

In the experiments, the local delamination phenomenon was observed as 
the predominant failure mode only for the n=4,6 and 8 specimens. The shear 
deformation model presented in this paper provides good agreement with the 
experimental data in this range. For n<4, edge delamination either in the 
mid plane or in the 25/90 interface was observed in the tests. Hence, the 
predictions of the local delamination models in this region are not of 
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consequence as long as they do not predict critical loads lower than those 
predicted by edge delamination models. Thus, it can be seen that the shear 
deformation model predicts the observed behavior with reasonable accuracy 
and can be used in conjunction with an appropriate edge delamination model 
to predict critical loads accurately for the complete range of n values. 
The edge de lamination model presented in References 6 and 12 can be used 
for this purpose. However, a separate model is required to account for the 
mid-plane (Mode I) edge delamination behavior. 

Conclusions 

A shear deformation model has been developed to analyze local delami- 
nations growing from transverse cracks in 90° plies located around the mid 
plane of symmetric laminates. The predictions of the model agree reason- 
ably with experimental data from [±25/90 n J s T300/934 Graphite Epoxy lami- 
nates. The predicted behavior is such that, in combination with an edge 
delamination model, the critical loads can be predicted accurately in the 
range of n from .5 to 8. 
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1 


Appendix A 

Sublaminate Analysis for Local Delaminations 

Interlaminar Stresses and Energy Release Rates 

A generic sublaminate is shown in figure 3 along with the notations and sign 
conventions. The interlaminar normal (peel) and shear stresses are denoted by P 
and T respectively with the t and b subscripts for the top and bottom surfaces 
respectively. The axial force resultant, shear force resultant and bending moment 
resultant are denoted by N, Q and M respectively. Plane strain conditions are 
assumed to prevail in the x — z plane and the thickness strain e zz is neglected. These 
assumptions lead to considerable simplification in the analysis. The displacements 
in the x and z directions are assumed to be of the form 

u = U(x) + z/3(x) (A.l) 

w = W(x) (A. 2) 

Here U represents the axial stretching and W is the transverse (thickness direction) 
displacement. This formulation recognizes shear deformation through the rotation 
(3. The equilibrium equations take the form 

= 0 


N, x + T t -T b 

Q,x + Pt ~ Pb 


o 


(A. 3) 
(A. 4) 


2 


M yX -Q + |(T, + T b ) = 0 (A. 5) 

where h is the thickness of the sublaminate. The constitutive equations in terms of 
the force and moment resultants are 


N — A n U tX + Bnfl tX 
Q = A 55 (ft + W <x ) 

M = BnU^ 4- D u ft yX 


(A.6) 

(A.7) 

(A.S) 


where A, B and D are the classical laminate theory axial, coupling and bending 
stiffnesses defined in the customary manner as 

An = J\C n dz 

h 

Bn = Cnzdz 


Dn = [\C n z'dz 
J 2 
r| 

A55 = J h C55 dz 

Here, the C s are the material moduli. For the case of plane strain in the x — z 
plane, the C s are defined as follows. 


(A.9) 


' 




f \ 

&XX 


C\\ Ci 3 0 


€-xx 

°zz 


Cl 3 C22 0 

< 

£zz > 

Txz 


1 

0 

0 

p 

Cn 


Txz 


The boundary quantities to be prescribed at the sublaminate edges are 

N or U 


M or ft 
Q or W 


3 

Farther, at the interfaces between sublaminates, appropriate reciprocal traction and 
displacement matching boundary conditions have to be used. 

The four sublaminates along with the loads acting on each axe shown in figure 4. 
Setting Pi and Ti as shown automatically satisfies the traction matching boundary 
condition at the 1-2 interface. From symmetry, we get w = 0 and zero shear stress 
along the bottom faces of sublaminates 2 and 3. This leads to w = 0 in sublaminates 
1,2 and 3. Thus, W has been prescribed in these sublaminates and the vertical shear 
force resultant Q cannot be prescribed at both ends of the sublaminates. Conse- 
quently, the calculated peel stress distribution will not be correct. In addition, at 
the 2-3 interface, the 0 s cannot be matched, since in these sublaminates, specifying 
0 is equivalent to specifying Q (through eq. A.7). Inspite of these simplifications, 
reliable energy release rate components can be estimated based on the interlaminar 
shear stress distributions. The mode I contribution can then be evaluated using the 
total energy release rate, which is not affected significantly by these simplifications. 

For the (±25/90 n ) s laminates under consideration, Bn is zero in all the four 
sublaminates. For sublaminates 1 and 2, the equilibrium equations and constitutive 
relationships can be written as 


£ 

H 

1 

II 

0 

(A.10) 

N 2tX + Ti = 

0 

(A. 11) 

Qi,x — P\ = 

0 

(A. 12) 

+ Pi ~ P 2 ~ 

0 

(A. 13) 

Mi, x + ^Ti-Qi = 

0 

(A. 14) 

M 2 , x + ^Ti-Q 2 = 

0 

(A. 15) 

H 

II 

f— * 


(A. 16) 




4 

N 2 — Au(2)C^2wc 

(A. 17) 

Q 1 = -^SS(l)/^l 

(A. IS) 

Qi = A 55 ( 2 ) 02 

(A.19) 

Mi = £„(,)&.* 

(A. 20) 

M 2 = D 11 (2)/?2 r r 

(A-21) 


The subscripts in brackets refer to the sublaminates to which the stiffness coefficients 
correspond. Equations A. 14, A. 15 and A. 12 can be rewritten in a modified form as 


M 1 , x + ^-N Ux = A 55(1) ft (A.22) 

M 2 ,X — ^N2,x = A 55 (2)/?2 (A. 23) 

Pi — Ql,x 

= M hxx + ^T ltX (A. 24) 

Matching the u displacement along the 1-2 interface implies 

Ul (“¥’*) = U2 (¥’ x ) 

or U t - = U 2 +!%/3 2 (A.25) 


Combining the equations to eliminate the displacement and interlaminar stress 
terms leads to the following homogeneous coupled system of ordinary differential 
equations. 





N 1 }X + N 2 ,x 

= 0 

(A. 26) 


Mi yXX + 

^Al,XX 

£>11(1) 

= 0 

(A. 27) 



t a ^- 

■^11(2) 

= 0 

(A.2S) 

iV, 

/z, M, 

n 2 

ho Mo 

= 0 

(A. 29) 

A n( i) 

2 -Dii(i) 

-4-11(2) 

2 £^11(2) 
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The solution is assumed of the form 


\ 

Ny 


' 4i ' 

n 2 

> = i 


Mi 


A 3 

m 2 

\ J 




>e' 


(A.30) 


Substitution of this solution into the governing equations results in the following 
system of algebraic equations. 

0 
0 




s 

0 


0 


0 

1 


L U(1) 




L n(2) 


2 _ ^ 55 ( 1 ) 
P H(1) 


0 

hi 1 


s 


2 D 


11 ( 1 ) 


2 _ 

•Dll(2) 

h"> 1 

TZ5 


11 ( 2 ) 



r 4i ' 


* \ 
0 

< 

4 2 

SX 

= 

0 


A 3 


0 




0 


(A.31) 


The corresponding eigenvalue problem has to be solved in order to obtain non trivial 
solutions. The eigenvalues turn out to be the roots of the following characteristic 
equation. 

s j-Bi-s 4 + B 2 s^ + B^ = 0 


(A.32) 


where 


By 

B 


1 

1 - 11 ( 2 ) 


1-11(1) 


1 

Di\(2) 


(t?) 

1 


2 = 


b 3 = 


1 


1 ..ffiglq. 

Ah( 2) -Dll(2) -^11(2) -^ll(l) 


4s5(l) , ^55(1) 

Unfit -Diini 77 


1 


^ll(l) 2- / ll(2) 
+ 


. - - ^5£ il + 1 - f > 55 l?J + j;55(2) 1 

-4-11(1) •C'll(l) 4-11(1) 1^11(2) B U (2) D n(i) 


§ 


1 4-55(1) 4-55(2) _ 1 4 55 (1) A 55 (2) 

4h( 2) -On(i) -Dll(2) 4-11(1) -Dll(l) £*11(2) 


For the material system and ply stacking sequence considered, B 2 > 4BiB 3 . Hence, 
the roots can be written as 


s = 0, ± 


— B 2 ± \J B 2 — 4 ByB 3 
2TT X 


(A. 33) 
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Only the zero and positive roots of eq. A. 33 are considered as they give exponentially 
decaying solutions, leading to finite values for the resultants at the sublaminate ends. 
Hence, the solution for N\ can be written as 


Ni = Oje* 1 * + a 2 e ,7X + oq 


Using this in eq. A.26 yields 


N 2 = —a ie 4ir — a 2 e'** + a 2 


(A. 34) 


(A.35) 


Substituting N x and N 2 in eqs. A. 27 and A. 28 provides the solutions for the bending 
moments as 


Mi = aikie SlX + a 2 k 2 e S7X 
hd 2 — aik 3 e six -)- o, 2 k^Q S3X 


The k parameters in the above solutions are defined as follows. 

¥4 


ki = 


'll(l) 


277 ^ _5 1 


k 2 = —Z- 


, 2 


fcs = 


TmT) $2 

¥4 

SS£21 _ 2 


^11(2) 
k < = -X 


¥4 


k'l 1 ( 2 ) 

If P is the applied force and w represents the specimen width, 


(A.36) 
(A. 37) 


(A.3S) 


(A. 39) 


(A. 40) 
(A. 41) 




(A. 42) 
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Using this in conjunction with eq. A.29 allows determination of the constants 
and a 2 . The following solutions for the stresses and the resultants can then be 
obtained. 


Nr 

_ .ill t - i P ^*ll(l) 

_ a t t + <.* +25 illl “+A 11(J) 

(A. 43) 

n 2 

- a r s »* n r a * x 1 P An(2) 

2u> A n (i) + A U ( 2 ) 

(A. 44) 

T\ 

= N hx 



= ajSje 4 * 1 + a 2 s 2 e SjI 

(A. 45) 

P i 

= M hxx + t$-T ltX 



= (fc + kt)ar$\e*' x + (k 2 + kj-)a 2 s\e* x 

(A. 46) 


The constitutive equations are used to write down the displacement solutions. The 
rigid body displacements of sublaminates 1 and 2 axe matched (in order to satisfy 
the displacement continuity condition) to obtain 


u ' ~ + & a 11(1) + (A - 47) 


u, = 


_ — r-^2 — e SJX + £■ 


L n(i)- S i 


•All(l) 5 2 


A 


11(1) 


+ ” 3 (A ' 4S) 


fix = + a 2 k 2 s 2 e a * x + + a 2 s 2 e s * x )} (A.49) 

fi 2 = a ^ — [ui^sSie* 11 + a 2 k4S 2 e 3iX + •^•(aiSie*' 1 -4* c^^e 521 )] (A. 50) 

•<*55(2) w 

The constants a x , a 2 and a 3 occurring in the solutions are determined using the 
boundary condtions. For sublaminate 3 the governing equations are 



= 0 

(A. 51) 

Qz,x + P 3 

= 0 

(A. 52) 

^3,x ~ Q3 

= 0 

(A. 53) 


£ 

cs 

II 

(A. 54) 

Ql — •d.55(2)^3 

(A. 55) 

= £>11 (2) /?3,x 

(A. 56) 


Matching U at the 2-3 interface and applying N 3 (a) = 0 gives 


N 3 = 0 


U 2 = U 2 (0) 


a i 


_£ 2 _ 


+ cl 3 


(A. 57) 
(A. 58) 
(A. 59) 


•51-^11(2) •S2-Ah(2) 

In order to solve for the bending moment, eqs. A. 53, A. 55 and A. 56 are combined 
to obtain 


M 3 ,xr - = 0 

■^ 11 ( 2 ) 

The solution of eq. A. 60 can be written as 

M 3 = <f> x sinhu; 3 :r -(- <j> 2 coshu> 3 x 
where the quantity u ; 3 is defined by 


(A. 60) 


(A.61) 


2 _ d. 55 (2) 
11 ( 2 ) 


^3 = XT 


(A. 62) 


(A. 63) 


Since the 0 matching conditon cannot be used at the 2-3 interface, the (remaining) 
boundary conditions are 

M 3 (a) = 0 
M 3 ( 0) = M 2 (0) 

The <f > s can be solved using the boundary conditions A. 63 as 

</f>2 = o, \k 3 + a 2 k^ 

h = 


— (j^cothu^a 


(A. 64) 
(A. 65) 
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The solution for sublaminate 3 can be completed by writing the following expres- 


sions. 


Qz = cosh o> 3 x -f <^ 2^3 sinhu^x 

(A. 66) 

fiz = [<^ 1^3 cosh u> 3 x + <^ 2^3 sinh u; 3 x] 

(A. 67) 

P 3 = tj— smhu; 3 x + <t> 2 coshu; 3 x] 
^1 1(2) 

(A. 68) 

The equilibrium equations for sublaminate 4 are 


£ 

H 

II 

o 

(A.69) 

Qa,x — 0 

(A. 70) 

M4 tX — Q\ = 0 

(A.71) 

The constitutive relations take the form 


N, — -4ll(l)C^4,r 

(A. 72) 

Qa — -4-55(l)(^4 + W 4tX ) 

(A. 73) 

Ma = -^ll(l)/?4,r 

(A. 74) 

Using eq. A.69 with the boundary condition N A (a) = yields 


II 

(A. 75) 

Similarly, using eq. A. 70 with Qi(a) = 0 results in 


<0 

Ox 

II 

0 

(A. 76) 

Matching M\ and M 4 at the 1-4 interface and using eq. A. 71 gives 


A^4 — d\}c j Cl 2^2 

(A. 77) 
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The U 4 displacement is obtained by integrating eq. A. 72 and using the displacement 
matching boundary condition U 4 ( 0 ) = £/i(0). 

u 4 = + + 33 ^) + “ 3 (A7S) 

Similarly, integrating eq. A. 74 and setting /? 4 (a) = 0 gives 


Pi = 7 ^ — [ai^i + a 2 k 2 )(x - a) (A.79) 

Using the solutions for Q 4 and /? 4 and the boundary condition W 4 (0) = 0 in eq. 
A. 73 yields the following solution for W A . 

W 4 = + a 2 k 2 }(^- - ax) (A. SO) 

In order to determine a l5 a 2 and 03, the following boundary conditions are used. 

JVi(0) = in 
MO) = MO) 

U\(—l + 0 .) = 0 


It is convenient to define the following parameters. 

*' = *!&<** + ¥> (*■“) 

*» - 1 (A - S2) 

* = 71^^+lt') ( A - S3 ) 

= rkt < A - S4 > 

&d = $3 — $1 ■+• ( 64 — 0 2 )a (A.S5) 

The nominal (far field) strain is given by 

„ _ P 1 


^ ^ 

Anjij 4- A U (2) 


(A.S6) 
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The a parameters are obtained as 

a i = ^ 11 ( 2 ) € \^ a(L (A. 87) 

a 2 = —^4.11(2) f (A.8S) 

u d 

«■ - < A - S9 > 

The specimen compliance C is defined as the ratio of specimen extension to applied 
load. This is obtained as 


r - 2U <( a ) 

= + “4 ( A£ >°) 

The total energy release rate associated with the crack (delamination) growth under 
a constant load P is given by 


Gf = 


_ P 2 dC 


tw da 


(A. 91) 


Using the compliance expression from eq. A. 90 in eq. A. 91 yields the following 
expression for Gj- 


r - P 2 

Gt ~s- 


(•An(i) in(i) + A U (2) 


(A. 92) 


where 


T — 1 -^- 11 ( 2 ) ^ 2^3 — f 1 — e _ 1 — e / A go) 

+ An(2) Aii(i) V s ! S 2 / 


T 1 ^n(2)(^3 + ^q)e- J>(< - a) --^ 1 +^Q)e- ,2(/ ~ a) 

2 ■4n( 1 ) + Au(2) 


(A. 94) 


The individual fracture mode contributions to the energy release rate can be cal- 
culated using the virtual crack closure method, based on the interlaminar stresses 
and displacements in the vicinity of the crack tip. From the assumed plane strain 
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condition, the mode III contribution is zero (G/j/ = 0). The mode II energy re- 
lease rate, G//, is calculated using the virtual crack closure technique while G/ is 
evaluated using 

G/ = G'/ — G// (A. 95) 

G// is calculated from the interlaminar shear stress and relative sliding displacement 
as 

G// = Hm J !Z\(z — 5)Au(i) dx (A. 96) 

In the absence of a singularity in the stress field, the result of the limiting process 
leads to the trivial result G// = 0. Hence, the limit is calculated as 6 tends to 
some finite value, say A. The value of A is chosen depending on the decay length 
associated with the problem i.e. the length within which the presence of the crack 
significantly alters the specimen response in comparison with the corresponding far 
field values. Evidently, the decay length in this problem is dependent on the eigen- 
values Sj and s?- The following value of A has been chosen in order to reasonably 
fulfil the decay length criterion. 



(A. 97) 


The relative sliding displacement An is based only on the difference U 4 — U 3 so that 
the kinematic condition of zero relative displacement at the crack tip is fulfilled. 
This also simplifies the calculations. The mode II energy release rate component is 
obtained as 

Gu = (A.9S) 

where J 3 and I 4 are defined as 


h = 


An(i) 


+ a ~ ) (ft ■ + §?) M 1 - - «-‘)] 


(A. 99) 
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Si A — 1 + e ,lA , s 9 A — 1 + e S2A \ + -^- 11 ( 2 ) 

a ' + ^ V in(„ 


(A. 100) 


Transverse Crack Spacing 


Shear Deformation Model 

The model presented so far has dealt with delaminations growing from a trans- 
verse crack. The same model can be modified to predict the spacing of these trans- 
verse cracks. In order to accomplish this, the delamination effect has to be isolated 
from the model. This can be achieved approximately by letting the crack length 
a tend to zero. This yields an approximation since the boundary conditions axe 
not accounted for properly by this limiting process. To get an accurate shear de- 
formation model, we consider only sublaminates 1 and 2 and apply the following 
boundary conditions for sublaminate 2. 

N 2 ( 0) = 0 (A. 101) 

M 2 (0) = 0 (A. 102) 

Using these boundary conditions in eqs. A. 37 and A. 44 yields two equations in a x 
and a 2 which can be solved to obtain 


_ k d P An( 2 ) 

1 fc 4 — k 3 2w An(i) + An(2) 

(A. 103) 

k 3 P Au(2) 

' x2 — k 4 2w An(i) + An( 2 ) 

(A. 104) 


The interlaminar shear stress can now be obtained using eq. A. 45. The saturation 
crack spacing corresponds to the distance from the crack where the broken plies 
regain their uniform stress/strain state i.e. where the interlaminar shear stress has 
decayed down to its far field (uniform) value. Practically, this distance is calculated 
by looking for the x where the interlaminar shear stress is some small fraction (say 
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.001) of its maximum value. The maximum shear stress evidently occurs at x = 0 


and is given by 


T i(mar) , 

1 — Gl$l + CL 2$7 


(A. 105) 


The crack spacing A can then be determined by solving the following transcendental 
equation. 


aT.Sie JlA + a?s?e ,2A 
GjSl + C2$2 


0.001 


(A. 106) 


Membrane Model 

A simpler model can be used to estimate the saturation spacing of the transverse 
cracks. This model treats the sublaminates as membranes i.e. the bending effects 
are ignored. The equilibrium equations for a generic membrane sublaminate are 


N s + Tt-T b = 0 (A. 107) 

%(T t -T b )-Q = 0 (A.10S) 

i 

The constitutive equations take the form 

N = ( A ' 109) 

Q = A ss /? (A. 110) 

The displacements are assumed to be of the following form. 

u = U(x) + z/3(x) (A. Ill) 

w = 0 (A. 112) 

The following governing equations can now be written 

Ni yX — Ti = 0 (A. 113) 

N 2 , i + T 1 = 0 (A. 114) 
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1 fT l -Q 1 = 0 

(A. 115) 

^Tx-g 2 = 0 

(A. 116) 

Nx = jiU\' X 

(A.117) 

Ni = 72^2,1 

(A.11S) 

Q 1 = ■^■55(1)01 

(A. 119) 

Q 2 = A 5 s(2)02 

(A. 120) 


(A. 121) 

where the 7s are defined as 


»2 

71 = Al, W uSS) 

(A. 122) 

72 = All(2 > iS 

(A. 123) 

Eqs. A. 113 and A. 115 can be combined as 


Qi = 

(A. 124) 

Using eqs. A. 119 and A. 117 in this leads to 


0i — _ 2 l ~a — ~l\Ui, xx 

Z -^55(1) 

(A. 125) 

Following a similar procedure for 02 yields 


& = hrr' u '** 

Z -^55(2) 

(A. 126) 

Using these two relations in eq. A. 121 leads to 


Mfc) Ar, u '~ = u '-&) tEsP'- 

(A. 127) 

Combining eqs. A. 113, A. 114, A.117 and A. 118 gives 


llUi iXX + 72 U 2 , XX = 0 

(A.12S) 
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Substituting this into eq. A. 127 results in 


U\' XX 






*■ 55 ( 1 ) \ * J ^- 55 ( 2 ) 

The characteristic roots of this differential equation are 


\yXXXT 4" i^Ul^xx 0 


s = 0, 0, ± 


n 


7i +72 


7i72 


3^ + (%) Z 1 


55 ( 2 ) 


(A. 129) 


(A. 130) 


The solution for U\ can then be written as 


U i — Aie six + A 2 X + Az (A. 131) 

where the As are arbitrary constants to be determined from the boundary condi- 
tions. The root S\ is the positive root such that a decaying solution is obtained in 
the negative x region. For the special case of 5 n (i) = Bn ( 2 ) = 0, the nonzero roots 
can be written in a simpler form as 

1 

— |_ -2 .. 

^ 55 ( 1 ) ^ 55 ( 2 ) 

The interlaminar shear stress can be obtained as follows. 


(A. 132) 


s 2 = jjylgl) 


Ti = N hx 

= 71 ^ 1 , XX 

= 7iAiSje iia: (A. 133) 

The maximum shear stress is 

Ti max) = ll A 1 sl (A. 134) 

Then, the saturation crack spacing A corresponds to 


e S| A = 0.001 


(A. 135) 
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Shear Lag Model 

This model allows for a nonlinear displacement field through the thickness of the 
sublaminate. Its fundamental assumption is that the shear deformation neglected 
in the classical theory of bending can be estimated using the shear stress. The 
sublaminate axial force equilibrium condition can be written as 

N x + (T t -T b ) = 0 (A. 136) 


The axial stress is assumed to be uniform and is given by 



The shear stress is estimated as follows 


(A. 137) 


&xz y z — 

' -A r .x 

= ~H~ 
= 

rt 


(A.13S) 


This can be integrated to obtain 

a xz = lL^ z + Ii^3± (A. 139) 

Neglecting transverse displacement, the axial displacement can be obtained by in- 
tegrating the shear strain, which in turn is obtained from the shear stress. 


u , = 


tT 

0-55 

1 


(Tt - T 6 )f + 21+^ 


u 


= U(x) + ^\(T t -T b )£ + (T t + T b )z 


(A. 140) 
(A. 141) 


where U(x ) is the mid-plane axial displacement. This displacement expression can 
be used to obtain an improved axial stress estimate as follows. 


® XT 


G 1 1 ^ ,X 


IS 


= c, 


n 


U' x 2CsS^ 1 ~ ^' x ~h "** ^ xZ 


(A. 142) 


The corresponding axial stress resultant can be written as 



(A. 143) 


The governing equations for the sublaminate are thus eqs. A. 136 (equilibrium), 
A*141 (displacement field) and A. 143 (constitutive relationship). Using these to 
model sublanunates 1 and 2 results in the following governing equations. 


Ni tX — T\ = 0 
-^2,x + T 2 = 0 

Ni = (7x1(1) 

N2 = Cn( 2 ) 


hiUi ' x - 

h? 


\h 2 U 2iX + ^cl (2 Ti, x 
Ui = Ul + 2 + 

U 2 = U2 + 2ck^[ Tl £ + TlZ . 

Displacement continuity at the 1-2 interface implies 

= u 2 (x,-^) 
or U 2 = U t - ^ 

Equation A. 146 can be rewritten as 

TT — Ni y hi m 

U '* - + 54C^ t -.' 

Combining eqs. A. 147, A. 151 and A. 152 results in 


AT (~* f h 2 N} h 2 Ty. x 

A2 " cn(2) l 


V 


hi j_ -J 12 . 


55(1) ^55(2) J 


(A. 144) 
(A. 145) 
(A. 146) 

(A. 147) 
(A.14S) 
(A. 149) 


(A. 150) 
(A. 151) 


(A. 152) 


(A. 153) 
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But from eqs. A. 144 and A. 145, we have 


N 2tX = -Ti = —N\ iX 


(A. 154) 


Using this in the differentiated form of eq. A. 153 leads to 


TuC 


2'~'11(2) 


+ 


ThC 


n(i)J 


Ni 


_ 1 


+ 


1 ^ 55 ( 1 ) ^ 55 ( 2 ) 


N- 


l, XXX 


(A. 155) 


The nonzero characteristic roots of this equation are given by 


2 _ o ( C ^ + ^2^11(2) A 

VAlCii(i)y V ^2^11(2) / V^2^55(l) + ^1^55(2 )/ 


(A. 156) 


This is the same as in the membrane model except for the factor 3 which is 4 in the 
membrane model. This difference is related to the fact that the axial displacement 
distribution through the thickness is parabolic in the shear lag model and linear in 
the membrane model. The crack spacing A for the shear lag model is determined 
as in the case of the membrane model but using the modified characteristic root. 



Table 1 Summary of Results 


number of 
90° plies 

Gt 

J/m 2 

Gii/Gj 

o- c 

MPa 

e c 

% 

1/2 . 

2.404 

0.276 

1313.9 

1.6747 

1 

6.752 

0.275 

784.0 

1.1685 

2 

22.849 

0.267 

426.2 

0.8058 

3 

51.049 

0.261 

| 

285.1 

0.6427 

4 

93.603 

0.256 

210.6 

0.5444 

6 

228.871 

0.250 

134.7 

0.4264 

8 

440.065 

0.247 

97.1 

0.3555 
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Fig. 2 Modelled Region and Sublaminate Scheme 
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Appendix III 

III.l Strain Energy Release Rate 

In this analysis, a delamination between belt and core sections is assumed to grow 
parallel to the belt direction in the tapered and uniform sections. These delamina- 
tions in each section are denoted by a and b respectively. The core section in the 
taper portion is modelled by two equivalent sublaminates. The stiffness properties 
are smeared to obtain the effective cracked and uncracked stiffnesses which are des- 
ignated by A u and A e as shown in Figure III.l. These stiffnesses change from one 
ply drop group to another with crack growth a by experiencing a sudden change 
at discrete locations. Therefore A u and A c can be represented in three consecutive 
regions as follows, 


• Region 1: 0 < a < l 


A» 


d + Zl — a 

d i / . / * / ~ a 

A-BD ^ A i A 2 A3 


(III.l) 


Ac = A 3 


• Region 2: l < a < 21 


A u 


d + 3/ — a 

d 1 l 2/ — a 

Ag D Ai A 7 


A <l b 

Ac = l Id THE 

A 2 ^ A 3 


(HI--) 


(HI. 3 ) 

( II 1-4 ) 


1 


o 


• Region 3: 21 < a < 31 


A u 

A c 


d ± 3/ — <x 

_il 4- -L 

Abu Aj 


a + b 


a-2t 

Ai 


+ 17 + 


l±b 

Ai 


(III-5) 

(III.6) 


where 

h = ply thickness 

d = length of uniform thick portion . 

/ = distance between two consecutive ply drop locations 

.4! = 6hQ 4 > - 2hQ° 

A 2 = 4 hQ 45 + 2 hQ° 

A z = 2h.Q 45 + 2hQ° 

A bd - 7h.Q° + 2hQ 45 
Q° = Qn of a 0 degree ply 
Q 45 = Q u of a ±45 degree ply 

Geometry of the sublaminate model is shown in Figure (III.l ) 
Also axial stiffnesses As, A a and Ap are given by 

d ± 3 / — cl 

A b = 


d , 3/— q 

Abd a bt 


(HI-7) 


3 



Figure III.l: Geometry of the Sublaminate Model 

Ap = .4 3 

= Abd 

w;here 

Abt =Taper belt stiffness 

For a membrane behavior, equilibrium equations are reduced to 

N, x = 0 

and the displacement field is assumed to be 

u(x,z) = U(x) 


(III.S) 
(III. 9} 


(III.10) 

(in.il) 


w = 0 


(III. 12) 
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The constitutive relations are represented by 


N = AnU* (III.13) 

The stress and displacement fields, are determined based on the stiffnesses derived 
in Equations(III.l-III.9). In this model, load is shared by the core and the belt 
portions according to their respective stiffness ratios 


Pi = 


PAb 

Ab + A u 


(III. 14) 


P 2 


PA, 

Ab + A„ 


( III. 15) 


where P is half of the total axial load applied at the ends. 

Using the Equations (III. 10), (III. 13), and the expressions for P] and P> from 
Equations (III. 14). (III. 15) the axial displacement at x = c can be written as 


PAbc P(d 4- 3/ -f b) / Ab Ab \ 
A,(A S + K) + (A b + Au) KAj, ~ TJ 

P(d + 3/ — a) / _ Ab \ 

{Ab A u) \ Abi / 


(III. 16) 


P A u c P(d + 3/ + b) / A u A u \ 
A 3 (Ab + A u ) (Ab + Au) \ A c Ap ) 

P{d + 3/ — a) / _ Au\ 

( Ab+Au ) V A c ) 


(III. 17) 


where Ab i is the belt stiffness in the pop-off region as shown in Figure III. 1 . 

A three-dimensional transformation is required in order to estimate the effective 
axial stiffness of the belt region Ab and Ap\- This is due to the belt layup and 
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the orientation of the different belt portions to the loading axis as shown in Figure 
III.l. The three-dimensional transformation is presented in section III.3. 

The tapered laminate is assumed to be fixed at x = 0. Therefore the external 
work done is given by 


W = PiU 5 + P 2 U e 


(III.18) 


Substitute from Equations (III.14) through (III.l 7) into Equation (III. 18) to get 


W 


P 2 (A B + A U ) 2 


(c — d — 31 — b) d" + (d + 3/ — a)(A B + A u ) 

(III.19) 


Al, A 2 s 


The strain energj r release rate G due to the external work done is determined 


bv 


G = 


1 dW 
2 P 2 dA 


(III. 20) 


where A is the delamination surface area. G is calculated for delamination lengths 
ranging from 0 to 60h. In the analysis, S2/SP250 Glass- Epoxy is used. Its properties 
are given in Table III.l. 


Table III.l: Material Properties of S2/SP250 Glass-Epoxy 


£n (MSI) 

E 22 (MSI) 

G 12 (MSI) 

G 13 (MSI) 

G 23 (MSI) 

IS 

7.3 

2.1 

0.87 

0.5 

] 

0.5 

0.275 
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III.2 Interlaminar Stresses 

In this part, an analysis for the interlaminar stresses in the belt-core interface in 
the tapered section will is developed. 

The simple analytical model assumes a beam model for the belt in the tapered 
section which is shown in Figure III.2 . Material and geometric discontinuities are 
modelled as extensional fc, : and concentrated shear springs < 7 ,- (£=1-4) as shown in 
Figure III.3. The resin pockets are assumed to be subjected primarily to shear stress 
and they are represented by a distributed shear spring with a constant stiffness G. 
The effect of the core is incorporated as elastic supports on the beam-belt model. 

A minimum complementary potential energy formulation is used to estimate 
the interlaminar stresses. The total complementary potential energy consists of 
bending, shear and extensional energy contributions, 


it = lb. + n, + n, + n. 


(111.21) 


where lit, II t . Il e , II*,. represent bending, shear and extensional energy components 
and energy stored in elastic springs, respectively. These are given as, 


n b 


1 r 31 M%s) 

2 Jo 1-^11 


(III. 22) 




*'=\i 


1 r 31 N 2 {s) 


ds 


>■11 


(III. 23) 


(III. 24) 


IU 


ir^i i3+ K + 

2 J 0 G, 2 k, 


I V . IL + IL + 1L + IL 

2k 2 2k 3 2 k A 2 gi 2 g 2 2 g 3 2 g A 


(III. 25) 




Figure III.2: Geometry of the Model 
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where Ri, T{ (t=l,2,3,4) are unknowns. The constant shear stress, c, due to resin 
filler is an additional unknown. The total number of unknowns in this formulation 
is nine. These unknowns axe constrained by following equilibrium equations. 

Ri = -Rz - 2 R a + 2 N u + N 22 - ^{Nn - N 2l ) (III.26) 

R 2 = 2R 3 + ZR 4 - 3 N u + Un u - N „ ) (III. 27) 

7\ = -T 2 -T z -T a - 3c/ - N x2 + N n (III.2S) 

where Nn, N 22 , N 2 1 and N 22 denote the components of the extensional load at two 
ends of the belt section. 

The bending moment, shear force and axial force in each of the three ply drop 
regions are written as 

• Region 1: 0 < s < l 

M(s) = —Ni 2 & + + R 4 S + Ta — (III. 29) 

V’(s) = Ni 2 — R 4 (III. 30) 


N(s) = N u -cs-T 4 


(III. 31 ) 
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• Region 2: l < s < 21 


M(s) = -N l2 s + + (R 3 + R<)s - R 3 l + ( T 3 + T 4 )| 


V(s) = N u -R 3 -R< 


N(s) = N u -cs-T 3 -T< 


• Region 3: 21 < s < Zl 


M(s) = -N 12 s + + (-R, + i? 3 + R a )s + (2R 2 - R s )/ + (T 2 + T 3 + T 4 )^ 


F(s) = iV 12 - R 2 -R 3 -R 4 


N{s) = N u -cs-T 2 -T 3 -T 4 


Therefore the bending energy in Equation (III.22) can be written as 


n fc = 


— f'\- 

Jo 


N X 2 $ + 4- R 4 S + T 4 - 


T2 


ds 


+ 2^ r i~ Nns + i s+ {r 3 + r * )s - r * 1 +(T 3 +T 4 ) t ~ 

1 [ 3l f . T ct 

+ 2dTJv {-^ 5+ 2 S + 


-R 3 - 2 R, + N u - -(N u - K, 


(111. 32) 

(111. 33) 

(111. 34) 

(111. 35) 

(111. 36) 

(111. 37) 


s 
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3i?3 + 6iLj — 6JV 12 + — (TV n — TV 21 )j l 

+(T 2 + T 3 + T 4 ) - | ds 


Similarly for the shear energy 


(III.38) 


n, — t (TV 12 — R 4) 2 ds + r £ (N 12 — R 2 — R 4 ) 2 0 

+ r £ l~ 2Nu + Rs + 2Ra + h {Nu ~ Nn ] 


ds (III. 39) 


where 


_ 3 


SGtA * 


The energy of extensional loads can be expressed by 


n e = 


1 r 1 


2A 


11 


I ( A n - « - T 4 ) 2 ds + A- f*(N n -cs-T 2 - J 4 ) 2 


ds 


+ T7~ A N " -cs-T 2 -T 3 - T 4 ) 2 ds 

2 An ^ 2 / 

The energy stored in the elastic springs is written as 


tt 3 c 2 , 1 

* 2 G 2 + 2 k 1 [ 


R 2 — 2R 4 + 2N U + N 22 — — (A'n — TV 21 ) 


+ 2^. 2 ^3 ^ 4 ~ } + "jM 2 — 3A r i2 + — (A T u — TV 21 ) 


(III.40) 
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Tl T 2 T 2 

2<?2 2^3 2<? 4 


(III.41) 


The complementary potential energy in Equations (III.38) through (III.41) is ex- 
pressed in terms of 6 unknowns, namely R 3 , R 4 , T { (i=2,3,4) and c. By minimizing 
these expressions the following system of linear equations is obtained 


/ 27t 2 l 3 Zl_ 9P\ 

v i 2 z> + g, + 77; 


tl 3 „ 5 tl 3 n 

C+ D Rz+ 2D R,+ 


8D 





T 3 


'9 t 2 l 2 

, 8 D 


5tl 3 

2D 


+ -)?<= + f t (Nn - JV») + ^(N u - N n ) 


t 2 P 

3 D 1 


(III. 42) 


tl 3 / 2/ 3 1 4 1 \ 

r^ C ' ! 'l~+7 _ 'j-T - -+- — I-R3-T 


D 


3 D fc 2 A* 3 


'«’-l + ±) Rt + UL T , 


,2 D A-j A 2 


4D' 


tl 2 


tl 2 


3 / 3 


■ , zZ—t — I — , 2 6 \ 2A r 22 

+ 2D T3 + 2 D Ta ~[2D + kr + k 2 ) Al2 + ~ 


+ 


( tl 3 t t \ , 

\12D 2k r l ~ k 2 l) ( Nn ~ N2 ^ 


(III. 43) 


5 tl 3 ( 3/ 3 2 6 \ /4Z 3 4 9 l \ tl 2 **/: 

2 B C+ \w + ^ + £) R * + ('d + ^ + £ + ^ r > 


3t/ 2 

2D 


^ ( 4l 3 4 9 \ Kr 2 N 22 tl 2 

T ' = \-D + k, + k 1 ) N " + — + ifl <">■ - "»> 


(III. 44) 
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'5 Pl 2 3 A tP 

JD + VJ C+ 4D 


R,+ 2 , £> iJ ' + (« + ^ + ^) 7i+ (^ + ^) T ’ 

/pi i \ tp x pi 

+ \Id + J 1 ) t ' = w n * + J, {Nn - N ^ + i d { n " ~ *”> (IIL45) 


W 3 A tl 2 n 5 tp (PI 1 \ /**/ J x 

, ^ + ^J C+ 2 J D il3+ 45'^ + (lD + ^) r2+ (^ + ^ + ^l r3 


4D 


/ Pi 1 \ 5<Z 2 . 1 Pl 

+ \2£> + 9l ) T * ~ I D Al2 + ^ + &D (iVl1 “ A ' 2l) (III. 46) 


9 PP zr 




ZtP 


PI 1 


( Pl 1 


8D 9l J - ' 2D * 2D— \4D 9l J - \2I> 5l / 

(ZPl 1 l\ 3^ 2 1 

pi 


+ 


16D 


(tfn-tf 2 i) 


(III. 47) 


Tlic concentrated norinal and shear forces at the ply drop regions and the inter- 
laminar shear in the resin filler are estimated by solving the simultaneous system of 
equations in (III.42) through (III.47) and using Equation (III.26) through (III.2S). 
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III.3 3-D Transformation of Stiffnesses 

It has been determined that a three dimensional transformation of stiffnesses is 
required in order to estimate the effective axial stiffness of the belt regions, Ag and 
Ag\‘ This is due to the belt layup and the orientation of the different belt portions 
to the loading axis as shown in Figure III.4. 

The loading axis corresponds to axis 1 in the 123 coordinate system which is the 
transformed system. The principal material coordinates are denoted by l',2’ and 
3'. 

The stress-strain relationships in the principal material coordinates for an or- 


thotropic laminate 

are given by 



•|y}6xl = [<?]sx6 {^xl 

(III.48) 

where 



<?n 

— (1 — ^23 ^32 ) 1 E'j\ 

(III. 49) 

Q 22 

= (1 - u zi u iz )VE 2 2 

(III. 50) 

Q 33 

— (1 ~ l/ 12 l/ 2l)l / E 33 

(III.51) 

Q 12 

= [y 21 + ^23^31 )!' E\\ — (^12 + ^13^32 -E"22 

(III.52) 

Q 13 

— (^31 + v 7 ,\ v 3 l)V E\\ = ( J/ 13 + ^23 U 1 2 )^ E 33 

(III.53) 

Q 23 

= { V Z2 + ^12 ^31 ) V E 22 = ( l/ 23 + ,y 21 ,y 13 E 33 

(III. 54) 

Q 44 

- G 23 

(III. 55) 

Q 55 

— G 31 

(III. 56) 

Q 66 

= Gi 2 

(III. 57) 
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Figure III.4: 

l = (1 - Z'i 2 J/ 2 l — Z' 2 3"32 — "31*'l3 — 2j' 12 /' 23 7/ 31 )' _1 (III. 5$) 

The presence of angle plies in the belt region making an angle 6 in the 12 -plane 


results in the following constitutive relationship 


M = [q] {*'} 

(III.59) 

where the transformed reduced stiffnesses Q tj are given in terms of reduced stiff- 
nesses Qij as 

Qn = c 4 Q n + 2 c 2 s 2 Qi2 + + 4c 2 s 2 Q 66 

(III. 60) 

Qn = s Qu + 2cVQ 12 + c a Q 2 2 + 4 c 2 s 2 Q $ 6 

(III. 61) 

Q\2 = c 2 s 2 Q u + (c 4 + s a )Q\2 + c 2 s 2 Q 22 — 4 c 2 s 2 Q 66 . 

(III. 62) 


Qee — ^ c ' s2 Qu ~ 8c 2 s 2 Q u + 4c 2 s 2 Q 2 2 + 4(c 2 — s 2 ) 2 Q &e 


(III. 63) 
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Q 33 

= Q 33 

(III. 64) 

Q 13 

= C 2 Ql3 + •S 2 Q23 

(III. 65) 

Q 23 

= -s 2 Qi3 + c 2 Q 2 z 

(III. 66) 

C 

= cos 6 


s 

= sinO 



Any ply in the belt portion of the taper makes an angle (3 with the loading axis 
if it is in the uncracked belt portion and an angle a if it is in the cracked belt 
portion. By performing a rotation about the 2-axis, the stiffness along the loading 
axis , takes the form 




W = \c\ { £ } 


(III. 67) 


where cr ^ and are in 123-axis system and C t7 represent the elements of trans- 
formed stiffness matrix in this coordinate system. 

Since we have assumed 


u(x,z) = U(x) 


(III. 68) 


and 


w = 0 


(III. 69) 


For plane stress condition in 1-3 plane (i.e. cr i2 = 0 ; i =1,2,3) stress strain 
relations reduce to 


^li — (C'xi — C12/C22) e n 


(III. 70 ) 



V 


17 

where 

C n = tQ u +2c 2 s 2 Q 13 + s < Q 33 + c 2 s 2 Q 5& (III.71) 

C n = c 2 Q 12 +s 2 Q 23 (III. 72) 

C 22 = Q 22 (III. 73) 

where c and s axe cosine and sine of the angle which the cracked and un cracked belt 
portions makes with the loading axis. 

The coefficient of e n in Equation (III.70) represents the transformed axial stiff- 
ness. This value is used in the derivation of Ab and Ab 1 . 


f 



